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ABSTRACT 


The problem of time difference of arrival (TDOA) is important in underwater 
acoustics for both passive and active sonar. Classical approaches to this problem are 
based on generalized cross-correlation (GCC) methods implemented in the frequency 
domain. After appropriate weighting of the cross spectral data in the frequency 
domain, an inverse discrete Fourier transform (IDFT) is performed and the peak of 
the resulting GCC function is located in the time domain. 

This thesis shows that the cross-spectrum of the data satisfies an appropriate 
signal subspace model; therefore the IDFT can be replaced with a signal subspace 
technique such as MUSIC. The result is an enhanced ability to locate the peak. Fur¬ 
ther, application of methods such as root-MUSIC or ESPRIT produce direct numeri¬ 
cal estimates for TDOA without the need to search for a peak. Results are presented 
for an extensive set of simulations using both synthetic signal data and data from 
a ocean acoustic propagation model (MMPE). Results are further presented for an 
application of the new method to target localization and tracking. In all cases results 
are compared using both the new methods and the classical methods. 
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EXECUTIVE SUMMARY 


In the passive sonar problem, signals received at two or more sensors or hy¬ 
drophones can be used to estimate the position and velocity of a detected acoustic 
source. Passive systems, unlike radar or active sonar systems, cannot control the 
amount of transmitted energy to be reflected from the source. However, the covert¬ 
ness of passive systems can be advantageous, both in military and biochemical appli¬ 
cations. In practice, the number of receiving sensors, the observation time and the 
ratio of the background noise to the source signal strength after propagation loss, 
when balanced against total system cost, dictate the feasibility of passive systems. 

In the ocean, sound usually arrives at each individual omnidirectional receiver 
through more than one path (e.g., direct and/or surface, bottom reflected paths). In 
order to deal with this problem from a signal processing point of view, it is useful to 
decouple two separate cases: the multipath and the planar problems. For multipath 
signals, a simple model of the received signal is that each receiver sees a signal plus an 
attenuated and delayed signal corrupted by additive uncorrelated noise, coming from 
different directions. For the planar problem (i.e., when all receivers and the source 
are in the same plane), it is usually assumed that the energy arrives at each receiver 
through only one propagation path in the same plane with all receivers and source. 
This means that the time delay to be estimated is the travel time of the acoustic 
wavefront between pairs of receivers, so that the source position and velocity can be 
estimated. 

In this thesis we deal with the time difference of arrival (TDOA) problem. 
We begin by examining previous research for estimating time difference of arrival 
known as generalized cross-correlation (GCC) methods. We refer to these as “classical 
methods.” With these methods, data from the two channels are transformed to 
the frequency domain to form the conjugate product XY* (cross-spectrum). After 
appropriate weighting, an inverse discrete Fourier transform (IDFT) is performed and 
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the peak of the resulting generalized cross-correlation function is located in the time 
domain. Our main goal in this thesis was to prove that for this kind of problem it 
is possible to use so-called signal subspace methods. More specifically, we showed 
that the aforementioned cross-spectrum of the data satisfies an appropriate signal 
subspace model. In this way we were able to replace the IDFT with a signal subspace 
technique such as MUSIC, Minimum Norm, PCLP, or similar procedure. In our 
problem we apply the subspace methods in the frequency domain and observe their 
results in the time domain. This is in contrast to their more usual application for the 
estimation of signals in noise. We also applied the subspace methods of root-MUSIC 
and ESPRIT, which produce direct numerical estimates for TDOA, without the need 
to search for a peak in the “pseudo-correlation” function. Both types of methods 
performed efficiently. 

We continued with a thorough series of simulations using synthetic data (MAT- 
LAB) and data created by an ocean acoustic propagation model (MMPE), each time 
testing the performance of all methods. The results showed that the subspace meth¬ 
ods performed as well or better than the classical methods, especially with low SNR 
conditions and with more difficult environmental data. Finally we applied the var¬ 
ious TDOA algorithms (classical and subspace methods) to the problem of target 
localization and tracking. All methods produced successful results. 
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I. 


INTRODUCTION 


The problem of time-delay estimation (TDE) is important in the field of under¬ 
water acoustics for both passive and active sonar. In this problem a signal is emitted 
from a source (e.g., a submarine) and is received at two spatially separated sensors. 
If si(f) represents the original undistorted source and ni(t) and n 2 (t) represent se¬ 
quences of uncorrelated, additive noise, then the signals received at two spatially 
separated sensors may be modeled as 

xi(it) = Si(t) + ni(t) 

X2 (t) = asi(t + D) + n 2 (t) (1.1) 

where D is the time delay between the two sensors. In a more complicated scenario, 
the received signals are attenuated and contain multipath reflections of the original 
source created by propagation through the ocean environment, in addition to the 
noise, as shown in Fig. 1. In the unknown source case (generally passive sonar), X\(t) 
and X 2 (i) and the ordinary cross correlation may be used to estimate the position 
and velocity of a moving source. For the known source case (generally active sonar) 
ordinary correlation involves using the original source and only one received signal 
to estimate the time it takes the signal to travel from the source to the receiver, i.e., 
m(t) = 0 in the model given above. In the absence of propagation distortion, this is 
matched filtering, which is the optimum linear method when the noise is white. 

Once the signal has been detected and the time delay D has been estimated, 
the time delay can be used to estimate the bearing angle B shown in Fig. 2. The 
bearing estimate is given by the approximate rule 

B “ cos ~\cD/L) (1.2) 

where c is the speed of sound in water, B is the bearing estimate and D is the time 
delay estimate. It can be shown that B is the angle that the hyperbolic “fine of 
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Figure 1. Model of direct and surface-reflected sound ray paths received at three 
sensors. 



Figure 2. Planar model of two receivers separated by a distance L. 

position makes with the axis of the receivers; hence the approximation Eq. (1.2) is 
increasingly accurate as the range to the acoustic source increases. 

The critical part for the passive bearing estimation problem, or in general 
for the localization, is the accurate estimation of time delay. In the literature, this 
problem is treated using terms like time delay estimation (TDE), time delay (TD), 
time difference of arrival (TDOA), group delay, time-of-arrival difference (TOAD), 
phase delay and others. For our purposes it is assumed that all these terms are 
identical and the terms TDE or TDOA will be used for the rest of our discussion. 
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A. PREVIOUS RELATED RESEARCH 


Two conceptual procedures are basically known according to [Ref. 1]. 

• An intuitive approach and a familiarity with detection theory. 

• A rigorous application of the maximum likelihood (ML) for white signals in 
white noise. 

In both conceptual procedures it is attempted to “advance” the delayed received sig¬ 
nal by a hypothesized amount in order to align it with the other received signal. In 
other words both received signals X\ (£), x-2 (£) contain the same transient S\(t), co ming 
from the same source, “buried” in different noise, rii(t), n 2 (£), and the only difference 
between them, is the time of arrival (TOA). After the above hypothesis either they 
are summed, squared and averaged as shown in Fig. 3 or multiplied and averaged 
as shown in Fig. 4. In both cases the hypothesized delays are adjusted in order to 



manmce J{0) 

Figure 3. Conceptual delay, sum, square and integrate configuration. 

maximize the configuration output. From the figures it is observed, that both config¬ 
uration outputs consist of “signal-cross-signal” terms. Further, both configurations 
are ML estimators for time delay under the assumptions that the signal and noises 
are white and mutually uncorrelated [Ref. 1]. When the signal and noise spectral 
characteristics are nonwhite, the received waveforms must be prefiltered with partic¬ 
ular equiphase filters (i.e., the prefilters must have the same phase characteristics) to 
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Figure 4. Conceptual cross correlator configuration. 


accentuate the frequency bands with good signal-to-noise ratio (SNR). It is of inter¬ 
est to note that a signal detector can be realized by comparing either configuration 
output to a threshold. Moreover in terms of detection (but not estimation) in the 
presence of noise, the system in Fig. 3 outperforms the system in Fig. 4 [Ref. 1]; 
however, the system in Fig. 3 requires prior knowledge of power levels in order to 
set the proper threshold. The system in Fig. 4 has a zero mean output in the signal 
absent, noise present case. 

A number of different approaches have been taken so that the conceptual 
systems in Figs. 3 and 4 can be achieved. The system in Fig. 3 can be implemented 
as shown in Fig. 5 which is called a time delay beamformer. It is presumed that s(t) is 
a sampled version of the original broadband time signal and that the sampling rate is 
large in comparison with the required Nyquist rate. In particular, the sampling rate 
is much greater than twice the bandwidth. The system configuration shown in Fig. 
4, including prefilters can be instrumented by a generalized cross-correlation (GCC) 
function, which in general is the inverse fast Fourier transform (FFT) of the product 
of a weighting function and the estimated complex cross-power spectrum . This means 
that data from the two sensors are transformed to the frequency domain to form the 
conjugate product XY*. After appropriate weighting, an inverse DFT is performed 
and the peak of the resulting generalized cross-correlation function is located in the 
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Figure 5. Time-domain beamformer implementation of one conceptual configuration 
for three hypothesized delays. 

time domain. All the methods that follow this general approach will be referred to in 
this thesis as classical methods and quite extensive research has been conducted on 
these methods [Ref. 2, 3, 4, 1, 5]. 

More recently, techniques based on higher order statistics such as bicorrelation 
and bispectra have been introduced to the problem [Ref. 6, 7, 8, 9]. These methods 
are claimed to have some advantage when there is Gaussian noise present which is 
correlated between sensors. In particular, since higher order cumulants of Gaussian 
noise are theoretically zero, these methods are practically “blind” to Gaussian noise. 

B. THESIS OBJECTIVES 

This thesis has the following goals. The first objective is to research and 
develop an algorithm that can be used for Time-Difference-Of-Arrival (TDOA) esti¬ 
mation based on subspace methods. The next objective is to evaluate the performance 
of this algorithm by comparing it with the respective algortihm based on generalized 
cross-correlation (GCC) methods implemented in the frequency domain, through a 
thorough set of simulations using synthetic (MATLAB) and model-based acoustic 
(MMPE) data. Finally the last goal in this thesis is to develop localization and 
tracking algorithms implementing the above methods, and to evaluate their perfor¬ 
mance in target’s extraction data (position-course-speed). 
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C. THESIS APPROACH 

In this thesis we suggest an alternative to the classical methods based on 
some more modern approaches to spectral analysis. In particular, we focus on meth¬ 
ods based on the signal subspace concept. We observe that the data product XY* 
computed in the classical methods satisfies a signal subspace model and therefore a 
number of well-known techniques such as MUSIC [Ref. 10], ESPRIT [Ref. 11] and 
others can be used to estimate the time delay. Mathematically, the procedure can 
be viewed as replacing the inverse DFT in the classical methods with one of these 
other techniques. As we later see from experimental results the new methods retain 
the accuracy of the classical methods, but give a much “cleaner” indication of the 
peak. Further, the use of methods such as root MUSIC or ESPRIT for the estimation 
provides direct numerical estimates for the time delay without the need to search for 
a maximum. We also have to mention that we have obtained similar behavior with 
the “Maximum Likelihood” or “Minimum Variance” method proposed by J. Capon 
[Ref. 12, 13]. Our focus in this thesis however is on the subspace methods. 

D. THESIS OUTLINE 

The remainder of this thesis is organized as follows. Chapter II provides a 
short introduction of the concepts of the classical methods and also describes the 
special characteristics of each one method, that is going to be used later in the sim¬ 
ulations from this family of methods. Chapter III discusses the general idea behind 
the subspace methods, gives a brief insight of all the methods from this group, and 
finally provides an argumentative presentation of the way that these methods are im¬ 
plemented into our problem (TDOA). Chapter IV presents the whole thesis approach 
to the TDOA problem, using various synthetic data (MATLAB software) and model 
based acoustic data using the Monterey-Miami Parabolic Equation (MMPE) model 
[Ref. 14, 15] for more realistic results, under different conditions/environments each 
time. Chapter V discusses the simulation results for both kinds of data (synthetic and 
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model-based acoustic) for the TDOA problem. Chapter VI develops a localization 
algorithm in two dimensions and discusses the implementation of the above methods 
in this localization problem in order to examine, how they ultimately function in the 
desired goal, namely the “Target Detection and Localization.” Chapter VII further 
extends the research to “Target Tracking Problem.” For this, one has to estimate not 
only target position, but also target course and speed, using Doppler measurements 
provided by MMPE. Chapter VIII presents conclusions and suggestions for future 
work. 
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II. CLASSICAL METHODS 


A. GENERAL APPROACH 

We have seen that the transmission of a short transient signal from a remote 
source monitored in the presence of noise at two spatially separated sensors can be 
mathematically modeled as 

*1 (0 = si(t) + ni(t) , 

x 2 (t) = asi(t + D) + n 2 (t ) , (II.X) 

where Si(t), ni(f), and n 2 (t) are real, jointly stationary random processes. Signal 
Si(t) is assumed to be uncorrelated with noise nj(t) and n 2 (t). 

One common method of determining the time delay D is to compute the cross 
correlation function 

Rx lX2 (t) = J r xx(t)x 2 (t - r)dt , (II.2) 

where T represents the observation interval. In order to improve the accuracy of the 
delay estimate D, it is desirable to prefilter X\ (t) and x 2 (t) prior to integration in Eq. 
(II.2). As shown in Fig. 6, X{ may be filtered through Hi to yield yi for i = 1,2. The 
resultant yi are multiplied, integrated, and squared for a range of time shifts, r, until 
the peak is obtained. The time shift causing the peak is an estimate of the true delay 
D. When the filters H\{f) = H 2 (f ) = 1,V/, the estimate D is simply the value of r 
at which the cross-correlation function peaks. 

The cross-correlation between xi(t) and x 2 (t) is related to the cross-power 
spectral density function by the Fourier transform relationship 

ik, ra (r) = /“ (II.3) 

J —oc 

When Xi(t) and x 2 (t) have been filtered as shown in Fig. 6, the cross-power spectrum 
between the filter outputs is given by 

<W/) = ffi U)HlU)G xm (f). (II.4) 
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Figure 6. Received waveforms filtered, delayed, multiplied, and integrated for a vari¬ 
ety of delays until peak output is obtained. 

where * denotes the complex conjugate. The generalized cross-correlation between 
£1 (t) and X 2 (t) is defined by 

w = /“ Mf)G nn (fy 2 ’ ,r df, (11.5) 

where 

Mf) = Mf)Hl(f) (II.6) 

and denotes the general frequency weighting. 

In practice, G XlX2 (f ) is not known a priori , and only an estimate G XlX2 (f ) of 

G x 2 x 2 (f) can be obtained from finite observations of xi(t) and x 2 (t). Consequently, 
the integral 

4m ( T ) = W)G X1X2 (f)^ fT df (II.7) 

is evaluated and used for estimating delay. Indeed, depending on the particular form 
of ipg(f) and the a priori information, it may also be necessary to estimate ipg{f') in 
Eqs. (II.5) and (II.6). For example, when the role of the prefilters is to accentuate 
the signal passed to the correlator at those frequencies at which the signal-to-noise 
(S/N) ratio is highest, then 4> g (f) can be expected to be a function of signal and noise 
spectra which must either be known a priori or estimated. 
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B. PROCESSOR WEIGHTING EVALUATION 


Before continuing to describe the weighting functions that are used in the 
generalized cross-correlation (GCC) in practice, it is informative to examine first the 
effect of processor weightings on the shape of Ry iy2 {T) under ideal conditions. For 
models of the form of (II. 1), the cross correlation of Xi(t) and a: 2 (f) is 

Rx iX2 (t) = aR SlSl (r - D) + Rn ltt2 (r). (II.8) 

The Fourier transform of (II.8) gives the cross-power spectrum 

<?—«(/) = + <W/). (II.9) 

If ni(t) and n 2 (t) are uncorrelated (G nin2 (f) = 0), the cross-power spectrum between 
Xi(t) and z 2 (f) is a scaled signal power spectrum times a complex exponential. Since 
multiplication in one domain is a convolution in the transformed domain, it follows 
for G nin2 (/) = 0 that 


Rx lX2 (r) = aR SlSl (t) © S(t - D ). (II. 10) 

where ©denotes convolution. 

One interpretation of (II. 10) is that the delta function has been spread or 
“smeared” by the Fourier transform of the signal spectrum. If Si(t) is a white noise 
source, then its Fourier transform is a delta function and no spreading takes place. An 
important property of autocorrelation functions is that |I? ss (r)| < R ss (0). Equality 
will hold for certain r for periodic functions. However, for most practical applications, 
equality does not hold for r ^ 0, and the true cross correlation (11.10) will peak at D 
regardless of whether or not it is spread out. The spreading simply acts to broaden 
the peak. For a single delay this may not be a serious problem. However, when the 
signal has multiple delays, the true cross correlation is given by 

Rnztir) = R SlSl (r) © £>£(* - A). (11.11) 
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In this case, the convolution with R SlSl (r) can spread one delta function into an¬ 
other, thereby making it impossible to distinguish peaks or delay times. Under ideal 
conditions where V/,(? Xll2 (/) = G XlX2 (/), / ip g (f) should be chosen to ensure a large 
sharp peak in Ry iy2 (j) rather than a broad one in order to ensure good t ime -delay 
resolution. However, sharp peaks are more sensitive to errors introduced by finite 
observation time, particularly in cases of low signal-to-noise S/N ratio. Thus, as with 
other spectral estimation problems, the choice of ipg(f) is a compromise between good 
resolution and stability. 

Having all the above in consideration, we are equipped with the appropriate 
background for the role that l ip g (f) is going to play. Thus, let us examine some 
generalizations of the cross-correlation function individually. 

1. Roth Processor 

The weighting proposed by Roth [Ref. 4], namely 1 

Mf) = G^XT) (IU2) 

yields 

Equation (11.13) estimates the impulse response of the optimum linear (Wiener-Hopf) 
filter 

*" (/) " CTI- 

which “best” approximates the mapping of « x (t) to x 2 (t). If m(t) ^ 0, as is generally 
the case for (II.l), then 


GxixAf) — G SlSl (f) + G nini (/), 


(11.15) 


and 


r(R) 

c yiV2 


(r) = S(r -D)® r - — ^df. 


Gsi.i (/) + Gnwi (/) 


(11.16) 


1 the subscript R is used here to distinguish the choice of ip g {f) 
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Therefore, except when G nini (f ) equals any constant (including zero) times G SlSl (f), 
the delta function will again be spread out. The Roth processor has the desirable 
effect of suppressing those frequency regions where G nini (f ) is large and G XlX2 (f) is 
more likely to be in error. 


2. Smoothed Coherence Transform (SCOT) 

Errors in G XlX2 (f ) may be due to frequency bands where G n2n2 (f) is large, as 
well as bands where G nini (f) is large. One is therefore uncertain whether to form 
ipR(f) = -pr—t T\ or iM/) = 77 — 777 ; hence, the SCOT [Ref. 3, 16] selects 


T X\X\ 


(/) 


Gx2X 2 (f) 

Mf) = 


yjG XlXl (f)G 

x 2 x2 (/) 

T his weighting gives the SCOT generalized cross-correlation 

J “OO 


(11.17) 


(11.18) 


where the coherence estimate is given by 


For H x {f) = 


sJgx^xJj) 


71112 (/) — 
and H 2 (f) = 


GxiXiif) 


\G XlXl 
1 


(H.19) 


\j G X2X2 (/) 


, the SCOT can be interpreted (see 


Fig. 6 ) as prewhitening filters followed by a cross correlation. When G XlXl (f ) = 
G X2X2 (f), the SCOT is equivalent to the Roth processor. If n\ 0 and n 2 7 ^ 0, the 
SCOT exhibits the same spreading as the Roth processor. This broadening persists 
because of an apparent inability to adequately prewhiten the cross power spectrum. 


3. Phase Transform (PHAT) 

To avoid the spreading evident above, the PHAT uses the weighting [Ref. 3] 


which yields 


V>p(/) = 


1 

KW/)r 



GxiXjjf) 

l<W/)l 


e’ 2nfr df. 


( 11 . 20 ) 


( 11 . 21 ) 
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For the model (II.l) with uncorrelated noise (i.e., G nxn2 (f) = 0), 


|GW/)| = aG^if). 


( 11 . 22 ) 


Ideally, when G XlX2 (f) = G XlX2 (f), 

Gxix 2 (f) jeu) _ pj2itfD 

KW/)| * 

has unit magnitude and 

flaW = «(r-D). 


(11.23) 


(11.24) 


The PH AT was developed purely as an ad hoc technique. Notice that for 
models of the form of (II.l) with uncorrelated noises, the PHAT (11.21) ideally does 
not suffer from the spreading that other processors do. In practice, however, when 
G X ix 2 (/) 7 ^ G XlX2 (f), then 6(f) ^ 2 tt fD and the estimate of R^ 2 (r) will not be a 
delta function. Another apparent defect of the PHAT is that it weights G XlX2 (f ) as 
the inverse of G SlSl (/). Thus errors are accentuated where signal power is smallest. In 
particular, if G XlX2 (f ) = 0 in some frequency band, then the phase 9(f) is undefined 
in that band and the estimate of the phase is erratic, being uniformly distributed in 
the interval [—tt, tt] radians. For models of the form of (II.l), this behavior suggests 
that ipp(J) should be additionally weighted to compensate for the presence or absence 
of signal power. The SCOT is one method of doing this. 


4. The Eckart Filter 

The Eckart filter derives its name from work in this area done in [Ref. 17]. 
Although it is not used in the experimental part of the thesis, derivations in [Ref. 
18, 19, 20] and [Ref. 21 ] are outlined here briefly for completeness. The Eckart filter 
maximizes the deflection criterion, i.e., the ratio of the change in mean correlator 
output due to signal present to the standard deviation of correlator output due to 
noise alone. For long averaging time T, the deflection has been shown [Ref. 19] to be 

d 2 = ^[/- 00 0 c^i(/)g2(/)g 3ia2 (/)d/] 2 m _ 

Ho l^i(/)l 2 l^(/)| 2 G nini (/)G n 2 n 2 (/)d/’ (IL25) 
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where L is a constant proportional to T, and G SlS2 (f ) is the cross-power spectrum 
between si(t) and s 2 (t). For the model (II. 1), G SlS2 (f) = aG SlSl (f)exp(j27rfD). 
Application of Schwartz’s inequality to (11.25) indicates that 

Hi(f)H’ 2 (f) = Mf)e+ 2 ’ r,D , (H.26) 


maximizes d? where 


i’E(f) 


aG 


SIS 1 


Gmn\ (. f)G 712712 (/)■ 


(11.27) 


Notice that the weighting (11.27), referred to as the Eckart filter, possesses some of 
the qualities of the SCOT. In particular, both weightings act to suppress frequency 
bands of high noise. Also note that the Eckart filter, unlike the PHAT, provides us 
with zero weight when G SlSl = 0. In practice, the Eckart filter requires knowledge 
or estimation of the signal and noise spectra. For (II. 1), when a = 1 this can be 
accomplished by letting 


Mf) = \G nn (f)\[G Im (f) - |GW/)|)' [G„*,(/) - l<W/)ll- 


(11.28) 


Since in our problem it was not possible to have knowledge or estimation of the 
signal and noise spectra, it was not possible to use this weighting processor for the 
experimental part of the thesis. All the above processors are summarized in Table I 
and can be justified on the basis of reasonable performance criteria, whether heuristic 
or mathematical. 


Processor Name 

Weight 

Mf) = u)mu) 

Cross Correlation 

1 

Roth Impulse Response 

1/Gx-iX^f) 

SCOT 

l/^Gx^ (f)Gx 2 x 2 (f) 

PHAT 

1/|Gzix 2 (/)| 

Eckart 

G s l s 1 f\Gn 1 7i 1 {f)G ri2 7 l2 (f)\ 


Table I. Candidate Processors 
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III. SUBSPACE METHODS 


A. CONCEPT AND PRINCIPLES 


An entire class of spectral estimates is based on the concept of signal and noise 
subspaces associated with the correlation matrix for a random process. Subspace 
methods apply primarily to locating discrete components (lines) of the spectrum. 
Pisarenko’s harmonic decomposition was the first of these methods and motivated 
other improved methods, such as MUSIC, that followed. These methods are all 
based on the fact that when the data consists of complex exponentials in noise, the 
frequency vector w, defined by 


1 

e ju 

e 5{N- 1)« 


(111.1) 


is an eigenvector of the correlation matrix. Its projection onto an orthogonal subspace 
complementary to the subspace defined by all of the signal vectors produces a null 
which can be exploited to estimate the frequency. The methods can be formulated 
either as a search for peaks in a function (called a pseudospectrum) or as a polynomial 
root-finding problem. In this section we give a short presentation of the general 
principles [Ref. 22] in these methods. 

Consider M independent signals in noise, the observation vector x, the noise 
vector rj, and the signals vectors Si, defined by N consecutive samples of the random 
process, and where it is assumed that M < N. Analytically, this means that the 
observation signals x[n] (transient & noise) are equal with 

M 

*[«] = £ S *N + n[n] , (III.2) 

t=l 

where the transients Sj[n] are given by 


Si[n] = A i e juJin . 


(HU) 
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Specifically, 


*[ 0 ] 
x[l] 

(III.4) 

x[iV - 1] 

#] 

. *?[!] 

V = | . , (III.5) 

vW - 1 ]. 

1 

e M 

(III.6) 

e j(N~l)oji 
A 

x = ^A i s i + r ? , (III.7) 

1=1 

where Si is defined by (III.6), and Ai is the complex amplitude of the i th signal, 

Ai = \M^. (III.8) 




This is the general signal model used in all of the subspace methods. If the noise is 
white, the correlation matrix is 


M 

Rx = 5Z PiSiSj T + all , (III.9) 

i=l 

where P { is defined by 

Pi = E{AiA*} = E{\Ai\ 2 } . (III. 10) 

The last two equations can also be written with more compact matrix notation as 

A\ 

a 2 

+ v (in.li) 

Am 
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and 


= SP 0 S* T + <r 2 I, 


(III. 12) 


where 


and 


S = 


Si S2 • • • Sm 



Po 


Pi 0 • • • 0 

0 P 2 ••• 0 


(III.13) 


(III. 14) 


[0 0 ••• P M \ 

It can be shown that the correlation matrix Rx has M eigenvectors lying in the signal 
subspace, which is spanned by the s i; and N — M eigenvectors lying in the orthogonal, 
complementary noise subspace. All N — M noise subspace eigenvectors correspond 
to eigenvalues A i = <j 2 0 while all of the signal subspace eigenvectors correspond to 
eigenvalues Aj > al . 

When the noise is not white, but is still uncorrelated with the signals, the 
foregoing development leads to the correlation matrix 


R* = SP 0 S* T + alV v 


(III. 15) 


instead of (III. 12). Here is a normalized covariance matrix that represents the 
covariance structure of the noise vector. In this case the whitening transformation 

y = S- x / 2 x (III.16) 


leads to the correlation matrix 


Ry = = TP 0 T* t + cr 2 1 , (III.17) 


where 

T = E“ 1//2 S (III. 18) 
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is the matrix with columns that are the transformed signal vectors 

t k = S- 1 / 2 s k ; k = 1,2,--- ,M . (III. 19) 

In the transformed vector space, the eigenvectors corresponding to the M 
largest eigenvalues span the signal subspace; those corresponding to the small est 
eigenvalues (equal to cr 2 ) span the noise subspace. The eigenvectors and eigenvalues 
satisfy the equation 

Ry e k = (S" 1 / 2 RxE“ 1/2 ) e k = A k e k (III.20) 

which can be written as the generalized eigenvalue problem 


Roc e k — AkS^e k , 


(III.21) 


where 

e k = S- 1 / 2 e k . (III.22) 

The basis vectors that span the signal and noise subspaces in the original coordinate 
system are 

b k = £j /2 e k = Ej /2 (Sj /2 e k ) (III.23) 

or 


b k = E„e k . (III.24) 

Neither the eigenvectors e k nor the basis vectors b k axe orthonormal in the usual 
sense. The eigenvalues satisfy the same conditions, however. Those corresponding 
to the noise subspace are the N — M smallest eigenvalues (equal to <r 2 ) and those 
corresponding to signal subspace are the M largest eigenvalues (all greater than cr 2 ). 

Let us now return to the case of white noise and the correlation matrix (III. 12). 
It is convenient for our later discussion to define the matrices of eigenvectors 


Esig = 


e i e 2 • • • eM 


(III.25) 
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and 


Enoise — 


®M+1 eM+2 


©N 


(III. 26) 


and also the two matrices of eigenvalues 

Xi 0 


A s ig — 


0 A 2 


0 0 


0 

0 

Xm 


(III.27) 


Am+i 

0 

... 0 


' ^0 0 

... o 

0 

Am+2 

... 0 


0 of 

... o 

0 

0 

l 

. * 


0 0 

• • • o-o . 


and 


Annlse — 


The complete matrices of eigenvectors and eigenvalues are thus given by 

E — [E sig E no i S e] 

and 

A = A “ z 

0 i 

Now observe that it is possible to write Rx as 

R x = EAE* t = E sig A sig E s *J + E noise A no i S e E * Q ise 
and to write R^ 1 as 

R ^ 1 = EA - 1 E* t = E^A^EJJ + E^A^JB;^. . 


(III.28) 


0 


(III.29) 


(III.30) 


(III.31) 


(III.32) 


Further, since the columns of E sig are orthonormal and define the signal subspace, 
this matrix can be used to form a projection matrix for the signal subspace, 

,-i. 


_ T7l ( T7V*T T7' \ T7l*T Tp T7»*T 

sig -^sig ^-^sig^sig J -^sig *^sig-*-' s ig > 


(III.33) 
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where the last step follows because E*jE sig = I MxM . Likewise, E noise can be used 
to form a projection matrix for the noise subspace 

Pnoise = E noise E*Q ise = I - P sig , (III.34) 

where the last equality follows because the subspaces are orthogonal complements of 
each other. 


B. TYPES OF SUBSPACE METHODS 
1. MUSIC 


The MUSIC (for Multiple Signal Classification) method forms a correlation 
matrix of some size N > M + 1 from which we try to find the eigenvalues and 
eigenvectors. Prom the previous presentation, the eigenvalues and the corresponding 
eigenvectors are divided into two categories; the ones that “belong” to signal subspace 
and the others to noise subspace according to the estimated eigenvalues. If the number 
of signals is not known, it can be estimated by looking at the smallest eigenvalues 
and finding the set that are approximately equal. This number is equal to N — M. 


The MUSIC method involves projection of the signal onto the entire noise subspace. 
For the MUSIC method we have the following approaches: 


1. Use the squared magnitude of the projection of w onto the noise subspace. 
Since each of the signals is orthogonal to the noise subspace, the quantity 
w Pnoisew = w* T E no j se E*T se w goes to zero for the values of the frequency 
where w = Si. The MUSIC pseudospectrum is defined as 


Pmu(^ w ) = 


1 


1 


W* T P no iseW 

and therefore exhibits sharp peaks at the signal frequencies where w = Si. 


w*T Enoise E£. seW 


(III.35) 


2. Use an alternative root-finding variation of the method called “roof MUSIC". 
In this approach an eigenfilter Ei(z) is defined as 

Ei(z) = ei[0] + eifl]*- 1 + • • • + e^N - , (III.36) 

where the e»[n] are components of the eigenvector ei. In this way the MUSIC 
pseudospectrum can be expressed as 


Pmu(^) = 


£,".m + i e<(z)e; (i/z*) 


(III.37) 


z=ei u 
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Since the denominator goes to zero at z = e ju)i (i = 1,2,- 
nator polynomial 

N 


Pm 1 u(z)= E BkMKWz") 


i—M+1 


• •, M), the denomi- 
(III.38) 


has M roots lying on the unit circle. These M roots (which are, in fact, double 
roots) correspond to the signal frequencies. 


A modification to the MUSIC method was proposed by Johnson and DeGraff 
[Ref. 23]. Prom the previous discussion, the MUSIC pseudospectrum can be 
written as 


Pmu(^ u ) = 


w 


•*T 


(E£m+ 1 e i e i T ) w 


(III.39) 


Since the eigenvalues for all the noise subspace eigenvectors are equal to a 
pseudospectrum for MUSIC that differs from (III.35) only by a constant can 
be defined as 


PmuW u ) ~ 


1 


w 


*T 


(ES: 


**T 


M+l e i e i ) w W 


*T 


(lZilM+1 A, e i e f T ) 


W 


(III.40) 


where the last equality follows because A, = a^,i — M + 1,- • ■ ,N. Equation 
(III.40) gives the pseudospectrum of the “modified MUSIC method 1” by John¬ 
son and DeGraff which differs from the first approach in that, in practice, the 
estimated eigenvalues are not all exactly equal. 


2. Minimum-Norm Procedure 

In the Minimum-norm procedure [Ref. 24, 25] we find a single appropriately 
chosen vector d in the noise subspace and define the pseudospectrum in terms of this 
vector as 

PMN(eP) = | w .T d |2 = w »T dd »T w * ( IIL41 ) 

The vector which lies in the subspace is chosen so that the squared magnitude ||d|| 2 
is minimized subject to the constraint that its first component is equal to 1. The 
resulting vector is given by 


d = 


-Er 


.rp ~-'noise* 
C* X C 


(III.42) 
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where Enoise is the matrix of noise subspace eigenvectors and c* T is the top row from 
the partition of E noise such that 


Enoise — 


,*T 


EL 


(III.43) 


If the components of d are denoted by d[0], d{ 1], 1], with d[0] = 1, then the 

rmnimum-norm pseudospectrum can be equivalently expressed as 


pMN{ei “ ] = \D$Pj \p ’ (IH.44) 

where D(z ) is the polynomial 

D ( z ) = E d[k]z~ k • (III.45) 

/:=0 

Thus the frequencies can also be found as the roots of D(z) on the unit circle. 

3. Principal Components Linear Prediction 

Another technique involving simple modifications to some of the previously 
presented methods suggests itself. Whenever there is added noise, let the estimated 
correlation matrix be represented in terms of only its eigenvectors and eigenvalues 
pertaining to the signal subspace. 

M 

Bi M) = EA,eier T . (III.46) 

i=l 

This is sometimes referred to as the principal components approximation of the core¬ 
lation matrix. This procedure, developed by Tufts and Kumaresan [Ref. 26, 27], 
amounts to eliminating some of the noise and increasing the overall signal-to-noise 
ratio. For model-based methods, a model derived from the reduced rank correlation 
matrix alone can then be generated and used to estimate the spectrum. 

Tufts and Kumaresan tried to eliminate the terms in the noise subspace using 
the pseudoinverse. This task could have even more successful results if a high-order 
correlation matrix was implemented. If the value of the prediction order variable P 
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was approximately in the same range with the available data samples ( N s ), then the 
rank of the estimated correlation matrix is automatically reduced. In order to have an 
exact equality between the rank of the estimated correlation matrix and the number 
of signals, then the following relationship should hold: P = N s — M/2, where P and 
N s are defined above and M is the rank of the pseudoinverse. This specific case is 
called Kumaresan-Prony case. It has the advantage of less computational demands 
since the SVD of the square matrix is the same as its eigenvalue decomposition, and 
there is no need for both estimations to be conducted. For minimizing the variance 
of the frequency estimates, however, a smaller value of P ~ |iV s is suggested [Ref. 
26, 27]. The filter coefficient vector is thus computed from 


m / '*T r \ 


(III.47) 


where R^ and r are the appropriate partitions of the correlation matrix Rx, and the 
ej, are the eigenvectors and eigenvalues of R x . The pseudospectrum is given by 


the final expression 


where 


PpcLp(e? w ) 


|w* T a| 2 


|A(e^)| 2 ’ 


(III.48) 


(III.49) 


If the n umb er of signals M is not known, it has to be estimated from the eigenvalues. 
The signal frequencies are then determined by either finding the peaks of (III.48) or 
by finding the roots of A(z). 


4. ESPRIT 

The ESPRIT (Estimation of Signal Parameters via Rotational Invariance 
Techniques) [Ref. 28, 29] method takes a somewhat different approach to frequency 
spectrum estimation by exploiting a certain invariance principle that exists for the 
subspaces of two overlapping data sets. Let us consider a data set of N + 1 samples 
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x[0], x[l], • • •, x[TV] and define the two vectors 



z[0] 


’ *[1] ' 

X = 

*[1] 

and x = 

x[2] 


T 

1 

H- 1 

1 _ 


x[TV] 


(III.50) 


Let B and B denote two matrices whose columns are basis vectors for the 
vector spaces associated with x and x. It can be shown [Ref. 22] that B and B' are 
related by a set of linear equations of the form 


B¥ = B' , (III.51) 

where & is an M x M square matrix and the eigenvalues of S' are of the form \ = e ? Ui , 
where cu* are the desired frequencies. 

In the next few fines there will be a short presentation of the steps in the ES¬ 
PRIT algorithm (TLS version) in order to have an understanding of its functionality. 
We start by defining the TV + 1-dimensional random vector x pertaining to TV + 1 
consecutive data samples x[0],x[l], • • • ,x[TV] and estimating the correlation matrix 
Rj from the data. Then we compute the eigenvectors and eigenvalues of R^,: 1 

R-x«fc = k — 1,2,--’ ,N + 1 . (III.52) 


If necessary, we estimate the number of signals M from the eigenvalues. Afterwards 
a basis spanning the signal subspace is generated and is partitioned as 


B = £rj 




(III.53) 


Since B and B' will not satisfy III.51 exactly, a Total Least Squares (TLS) approach 
to computing’T is used. The details of TLS are beyond the scope of this presentation. 

x If the noise is not white III.52 can be replaced by a generalized eigenvalue problem. 
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However, a general outline would be to first compute the matrix V of right singular 
vectors of 

B B' 

and then partition V into four M x M submatrices 


(III.54) 


V = 


Vu V 12 
V 21 v 22 


The desired estimate for T' is then given by 


(III.55) 


*TLS = -y 12^22 • 


(III.56) 


Finally, the desired frequencies are computed as the angle of the complex eigenvalues 
of T'xls as shown below 


u k = /At; k = 1,2,---,M. (III.57) 

C. APPLICATION TO THE TDOA PROBLEM 

A model for the TDOA problem as we have already seen is as follows. We 
assume that a short transient signal d[n) is emitted from the source. We further 
assume the signal received at each sensor is subject to amplitude attenuation and 
additive noise so that the signals x and y can be written in discrete time as 

x[n] = d[n] + io[n] 

y[n] = Ad[n — L]+u[n\, (III.58) 

where A is the relative amplitude, L is the time delay to be estimated, and u and w 
are the additive noise terms. Note that we assume the transient signals, d, received 
at each location are identical but arrive at different times. This assumption will be 
relaxed in the simulations presented later [Ref. 30]. Let us define R[k] as the product 

R[k] = X[k]Y*[k), (III.59) 
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where -U[/c] and Y[k] are the DFT’s of the sequences x[n] and y\p\ computed at a 
size equal to at least twice the length of the data. Thus the inverse transform of R[k] 
represents an estimate of cross-correlation in the time domain. Now let u 0 represent 
the frequency resolution of the DFT. The DFT sequences X[k] and Y[k] are then 
given by 

X[k} = D[k] + W[k] 

Y[k ] = AD[k]e~ iuJokL + U[k] (III.60) 

and from (III.59) and (III.60) 

R[k] = A\D[k]\ 2 e? u ° kL + {D{k]U*[k} + AD*[k)W[k)^° kL + W[k)U*[k}}. (III.61) 

Let us also assume that the signal d[n] is broadband, and that the source spectrum 
is flat and therefore the magnitude |D[/c]| is approximately constant. Further, if 
the observed time sequences (x and y) are sufficiently long, the points in the DFT 
sequence for the term enclosed in brackets will be approximately uncorrelated. Thus 
the sequence satisfies the basic model for signal subspace analysis [Ref. 22], but 
with time and frequency interchanged. 

For the following discussion, we shall consider both related approaches to 
estimating TDOA. In the first approach a steering vector of the form 

|"l g) u 0 l g)2 u 0 l ' ' _ gj(Ar-l)a; 0 2j 

is projected onto the noise subspace and the reciprocal of the result is plotted as 
the lag parameter l is varied. (Here N is the dimension of the observation used in 
the chosen subspace method.) The resulting function, which we refer to as a delay 
indicator function (DIF), peaks at the desired estimate l = L. Although the DIF 
resembles a generalized cross-correlation function, it is strictly speaking not a GCC. 
It is analogous to the pseudospectrum used in the spectral estimation problem. In 
the second approach, direct numerical estimates are computed and no DIF is formed. 
Root MUSIC and ESPRIT follow this approach. 
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IV. THESIS APPROACH TO THE 

PROBLEM 


In order to verify all the theoretical results described in Chapters 2 and 3, var¬ 
ious simulations were conducted. All the simulations are divided into two categories 
depending on the kind of the data that were used. The first category used synthetic 
data generated from equations in MATLAB while the second category used data gen¬ 
erated by the Monterey-Miami Parabolic Equation (MMPE) acoustic propagation 
model [Ref. 14, 15]. 

All the simulations using the synthetic data assumed the same shape and form 
and length of transient at each receiver but for each one we added different noise with 
variations according to kind (white-colored), characteristics (variance), and amount of 
noise (SNR). White or colored Gaussian noise was generated in MATLAB and added 
to the signals to achieve a specified signal-to-noise ratio (SNR). The SNR specifically 
was defined by the formula 

_L v^-i s 2r„i 

SNR = — Li ; (IV. 1) 

cr N 

where N is the number of the samples of the signal, s[n], and <j 2 n is the variance of 
the additive noise (white or colored). For the MMPE data, the received transients 
were different, in general, and noise was added with variations as before. In this way 
we kept all the basic assumptions that we made in Chapter 1 for the TDOA problem. 


A. SYNTHETIC DATA 

For the synthetic data five different kinds of transients were used: 

1. exponential transient of the form 

$ u\ _ ( exp(-df), 0 < t < T 
^ ' I 0, otherwise ’ 


(IV.2) 


where d is the decay factor and T is the time duration (length) of the expo¬ 
nential transient; 
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2. sinusoidal transient of the form 


s(t) 


= { si 


sin(27 vft), 0 < t < T 
0, otherwise ’ 


(IV .3) 


where / is the frequency and T is the time duration (length) of the sinusoidal 
transient; 

3. damped sinusoidal transient which is a combination of the first two transients, 


s(t) = 


exp(-dt) • sin(27r/t), 0 < t < T 
0, otherwise 


(IV.4) 


4. linear swept-frequeney cosine generator ( chirp(T , fo,T u f x )) transient, 


_ | e ( j2n { fot+ i 0t2 }), 0 < t < T 
( 0, otherwise 


(IV. 5) 


This function creates samples of a linear swept-frequeney cosine signal at the 
time instances defined in array T, /o is the instantaneous frequency at time 0, 
and fi is the instantaneous frequency at time T x . Both f 0 and fi are in Hertz. 
In our case the instantaneous frequency sweep fi{t) given by 


where 



fo "i - 0 < t <• T 

0, otherwise ’ 


0 = 


fi ~ fo 
T x ’ 


(IV. 6) 


(IV.7) 


5. 


damped chirp transient which is the combination of the first and fourth tran¬ 
sients, 


s 



exp(-dt) • e( i27r [ /ot+ ^ t2 ]), 
0 , 


0 < t < T 
otherwise 


(IV.8) 


In this way we were able to compare the performance of the classical and 
subspace methods not only by varying the characteristics of the-noise but also by 
changing the characteristics of the original transients. This was done in order to have 
a more complete evaluation of the above methods. 
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B. MODEL BASED ACOUSTIC DATA 


Although they provide a way to test algorithms under ideal conditions, we 
understand the simulated data in MATLAB does not approach the real conditions 
that occur in the ocean. For this reason we used data from the MMPE model to 
test the TDOA methods under more realistic conditions. This allowed us to test 
them in more complex environments, where additional factors like area geography, 
water column and bottom characteristics, sound speed profile, etc., can change and 
influence the result of the sound wave propagation. 

1. Monterey-Miami Parabolic Equation (MMPE) Model 
and Split-Step Fourier (PE/SSF) Algorithm 

Before presenting data generated by the MMPE model, it is informative to 
describe what the model is what the model is and what it is capable of. We begin by 
defining the time-harmonic acoustic field in cylindrical coordinates (r, 2 , 4>) 

P(r, 2 , <f>,ut) = p(r, z, 0)e -lwt . (IV.9) 


The result of the combination of Eq. (IV.9) with the wave equation, also in cylindrical 
coordinates, produces the Helmholtz equation, 


IJL ( r ?i) + 

r dr l dr ) r 2 d(f) 2 dz 2 


1 9 P + 7 T? + k 2 n 2 (r, 2 , <f)p = -4nP 0 8(x - x " s ) , 


(IV. 10) 


where k 0 = ^ is the reference wavenumber, n(r, 2 , <fi) = ^ is the acoustic index of 

refraction, c 0 is the reference sound speed, and c(r, 2 , <j>) is the acoustic sound speed. 
Most of the properties of the environment are provided by c(r, 2 , 4>). Regarding the 
source function, it is assumed to be a point source located at a distance of r = 0m 
and depth z = zs and with reference source level P 0 . This reference source level is 
the pressure amplitude at a reference distance of R 0 = 1 m, and 


J((x)) = 


1 

27r r 


8(z- z s )8(r) 


(IV.ll) 


is the Dirac-delta function defining the point source contribution. 
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For simplification of the Helmholtz equation, we assume that the cylindrical 
spreading dominates the propagation, and the magnitude of the pressure vector will 
be given by 

p(r, z ) = -j=u(r,z). (IV. 12) 

Substituting Eq. (IV. 12) into the Helmholtz equation and neglecting the source term, 
we end up with the following equation 


d 2 u 1 d 2 u d 2 u 1 2 ( 2 i 1 

dr 2 + ~?d^ + di 2 + M n + 4fcP 


u = 0. 


(IV. 13) 


In Eq. (IV. 13), we observe that both the final term and the second term , which intro¬ 
duces azimuthal coupling between different radials, are proportional to —. Therefore 

r 2 ! 

invoking the uncoupled azimuth approximation, these terms will be neglected in this 
analysis. 

Defining the operator notation 


P — — 
op dr 


(IV. 14) 


where 


Qop — (// + £ + l ) 1 ^ 2 , 


2 , Id 2 

£ = n - 1 ’ ^Md? 


(IV. 15) 


(IV. 16) 


the homogeneous form of Eq. (IV. 10) then becomes 

(P2 p + k 2 o Q 2 op )u = 0 . (IV. 17) 

Proper factorization of the outward propagation field is obtained by defining 

u = Q- 1/2 *. (IV. 18) 

Because of the small range dependence of the environment it is also assumed that the 
commutator [Pop, Qop] is negligible. The outgoing wave then satisfies [Ref. 14, 15]. 


- i K 1 ^; = Qo P *. 


(IV. 19) 
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The physical meaning of the Eq. (IV. 19) is the representation acoustic energy 
in the waveguide, particularly the forward propagating acoustic energy when the 
backscattered energy is not significant. In other words Eq. (IV.19) is considered 
to be the basic equation upon which all one-way PE models are “built”. The next 
step is to create a method for generating solutions to this equation by developing 
approximations to the pseudo-differential operator Qop. 

For the creation of the numerical algorithm which will solve the PE, we observe 
that the acoustic field can be divided into a slowly modulating envelope function 
and a phase term which oscillates at the acoustic frequency. The PE field function, 
ip(r, 2 , (/>), is defined as 

= ipe ik ° r (IV.20) 

or, in terms of the acoustic pressure, 

p(r, 2> 0) = z. <t>)e ikor . (IV.21) 

v T 

The last equation is scaled in such a way that at r = R 0 , \ip\ = 1 and |p| = P 0 - Com¬ 
bining Eq. (IV.21) with the Helmholtz equation will provide us with wave equation 
for the PE field function, 

^ = -ik Q ip + ikoQopi 1 = -ik 0 Hopi>, (IV.22) 

or 

where 

Hop = 1 - Qop (IV.23) 

is a H amil tonian-like operator which defines the evolution of the PE field function in 
range. 

In Eq. (IV.22), the function ^ is a vector (in z) in Hilbert space. The rela¬ 
tionship between the values of ip at different ranges can now be expressed as 

ip(r + Ar ) = $(r)ip(r). (IV.24) 
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To propagate the solution out in range requires a representation of the propagator 
$(r). The MMPE model uses the split-step Fourier (PE/SSF) method to compute 
PE solutions. 

In the SSF algorithm the operators Hop and are a combination of scalar 
and differential operators. For the appropriate functionality of the SSF algorithm 
the different terms within H ^ should be separated, requring an approximation to the 
square-root operator. In the MMPE model, the wide-angle PE (WAPE) approxima¬ 
tion of Thompson and Chapman is employed [Ref. 14, 15], such that 


Hop ~ Top + Uop 


(IV. 25) 


where 


and 


Top = 1 — 


1 + 


1 

k%dz 2 


1/2 


(IV. 26) 


Uop = ~(n - 1). (IV.27) 

In this form, the differential operator has been separated from the index of refraction 
term as required for implementation with the SSF technique. In 2 -space, the operator 
Uop is a simple scalar multiplication operator, in other words a diagonal matrix, 
but the operator Top is not a diagonal matrix, so different depth eigenfunctions are 
coupled. However, in vertical wavenumber space, the corresponding operator Top is 
diagonal. 

The MMPE propagator function is then based on the “centered step” scheme, 


<£( r ) = e -iko^U op (r+Ar) e -ik 0 ArT 0 p e -ik 0 ^U op {r) 


(IV.28) 


where error analysis shows that this scheme provides third order accuracy in A r, and 
is the method used in the MMPE implementation. The general algorithm behind the 
PE/SSF implementation is then as follows. The PE field function ip is specified at 
some range r in the 2 -domain. A multiplication of the 2 -space operator e~ iko ^ u °p^ 
defined at the beginning of the range-step is applied. A transformation is then made 
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to the domain followed by a multiplication of the -space operator e - lk ° ArT °p . The 
result is then transformed again to the 2 -domain followed by a multiplication of the 
z-space operator e~ lko ^ :u °p( r + Ar ) defined at the end of the range-step. The final result 
is the field function at r + Ar. The discrete fast Fourier transform (FFT) subroutine 
employed in the numerical code assumes the convention 

1>(z) = FFT(ip(k z )) (IV.29) 


and 

iKfe) = IFFTMz)). (IV. 30) 


Therefore, the PE/SSF implementation can be represented by 

ip(r+Ar,z) = e - ik o^Uo P (r+Ar,z) x FFT{e- ik ° ArT °p (k * ) x I FFT [e _ifc °^^ (r ’ z) x ip(r, z )]} , 

(IV.31) 


where, in fc z -space, 




(IV.32) 


2. Data Modulation-Demodulation 

The MMPE model is a broadband, full-wave acoustic propagation model based 
on the parabolic approximation to the Helmholtz equation. After defining all the nec¬ 
essary parameters for the entire environment (such as sound speed profile, bathymetry 
of the water/bottom interface, the acoustic parameters of the bottom, the source in¬ 
formation), this model computes the complex pressure, transmission loss, and arrival 
time structure for any point on a computational grid in range and depth. 

Broadband results are obtained by running the MMPE model for all discrete 
frequencies in the chosen bandwidth. The travel time results are then realized by per¬ 
forming a Fourier synthesis and a necessary multiplication by some window function, 
S(f). The complex arrival structure of the pressure field can then be written 

/ oo P D roc 

S(/)p(r, z, f)e~ i2nft df = -*£■[ 

(IV.33) 
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By defining the reduced time T = t - (^), the phase factor e~ ik ° r = e~ i2nf t can be 
absorbed, such that 

P P r°° 

P(r, z, T ) = J ^ S(f)ty(r, z, 0)e~ i27r/r d/. (IV.34) 

Note that the use of reduced time does not influence the autocorrelation function. 

For all the environments tested the same window function was used. The 
modeled transmitted pulse had a center frequency of 400 Hz with a bandwidth of 
128 Hz and with 256 frequency bins computed in the DFT. This creates a signal 
of duration 2 sec in the time-domain. In order to obtain the real pressure for each 
sonobuoy (receiver) while maintaining the initial assumption that all the receivers 
were receiving the signals of the same duration, we developed the following procedure. 

• First the frequency-domain complex pressure vector, corresponding to each 
receiver, was multiplied by a “Hanwin” factor. This “Hanwin” factor is a 
Harming window with length equal of the number of the frequency bins. 

• In order to add the carrier in the baseband signal, we had to append to both 
sides of the frequency-domain vector a number of zeros, creating a total length 
of 1024 bins. This will assist us with the frequency resolution that we are going 
to work with in order not to lose any information from the original transient. 
In this way, that part of spectrum corresponding to positive frequencies is 
formed. 

• The part of the spectrum corresponding to negative frequencies was formed 
by taking the conjugate reverse of the original vector and appending it to the 
original vector. This produced a complete frequency-domain representation of 
the measured signal over 2048 frequency bins. The corresponding sampling 
frequency on the receiver was then 512 Hz. 

• Performing an FFT on this complex frequency-domain signal produced a real 
time-domain signal that represented the measurement on the receiver. Due to 
a small imaginary component due to the numerical processing, the real part 
of this signal was extracted for use in the TDOA algorithms. 

• Before adding noise to the transients, all transients must be in absolute time 
and not in reduced time, as in the previous step. For this procedure, a number 
of zeros are pre-appended, which for each pair of receivers is equal to the 
fraction of the difference of the two distances (of the 2 receivers) divided by 
the factor dt = 1/1024. 
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3. MMPE Environments 

Before proceeding with the simulation results, let us describe the four test 
cases that were used with the MMPE model. This will provide a common reference 
not only for the TDOA simulation results but also for the localization and tracking 
problems presented in the next chapters. The four environments are: 

1. Flat Bottom 

2. Sound Channel 

3. Shelf Break 

4. Internal Waves 

The first two of these are range-independent (RI) while the last two are range- 
dependent (RD). The main difference between these two categories is that in the 
range-independent case the received signal at each sonobuoy depends only on the 
distance between the source and the receiver. In other words, the sound speed profile 
and bottom are the same at every point in the area of interest. In the range-dependent 
case the sound speed profile and/or bottom changes according to range and direction 
from the source. This makes it more difficult to simulate since one has to be aware in 
more detail about the environmental conditions. Range-dependent simulations pro¬ 
vide more realistic data consistent with the complicated conditions that occur in the 
ocean. 

a. Flat Bottom 

This environment is defined by a flat ocean bottom and isospeed water 
column. In our case we have regularly spaced values of bottom parameters. At each 
range location, the bottom is assumed to be a homogeneous half-space below that 
position. The specific values of the properties for the case are listed in Table II: 

b. Sound Channel 

This range-independent case is similar to the previous one with the 
following major differences. First, the bottom bathymetry is varying and the sound 
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Bathymetry 

flat (300 m depth) 

Sound Speed Profiles 

isospeed (1500 m/s) 

Water Density 

1 g/cm 3 

Water Attenuation 

0.0 dB/A 1 

Bottom Sound-Speed 

1700 m/s 

Bottom Density 

1.6 g/cm 3 

Compressional Bottom 
Attenuation 

0.1 dB/km/Hz 

Shear 

None 


Table II. Environmental Properties for RI Test case: “Flat Bottom”. 


speed profile, which is the same for the entire computational grid, is defined by the 
following equations: 


c= < 


1515 4- 0.016z 


1490 


1 + 0.25 50°” ) q. 


when z £ \z sur f , Zd uc t\ 

r (IV. 35) 

when Z £ induct , ^max] , 


where z aur f — 0m is the surface depth, z<i UCt — 75m is the depth of the surface duct, 
Zmax = 500m is the maximum depth and z axis = 200m is the depth of the axis of the 
sound channel. The sound speed profile is shown in Fig. 7. The specific values of the 
properties for the case are listed below in Table III. 


Bathymetry 

varying 

Sound Speed Profiles 

varying 

Water Density 

1 g/cm 3 

Water Attenuation 

0.0 dB/A 

Bottom Sound-Speed 

1700 m/s 

Bottom Density 

1.6 g/cm 3 

Compressional Bottom 
Attenuation 

0.1 dB/km/Hz 

Shear 

None 


Table III. Environmental Properties for RI Test case: “Sound C hann el” 
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Figure 7. Sound Speed Profile for RI Test case: Sound Channel, 
c. Shelf Break 

The next case has to do with a range-dependent environment similar 
to that found in the regions of continental shelf breaks. This environment is defined 
by both a varying water column sound speed and bottom bathymetry with homoge¬ 
neous bottom acoustical properties. This environment is loosely based on shelf break 
front structures observed during the ONR Primer experiment [Ref. 31]. The general 
properties of this test case are defined in Table IV. 


Bathymetry 

varying 

Sound Speed Profiles 

varying 

Water Density 

1 g/crfr* 

Water Attenuation 

0.0 dB/A 

Bottom Sound-Speed 

1700 m/s 

Bottom Density 

1.5 g/cm 3 

Compressional Bottom 
Attenuation 

0.1 dB/km/Hz 

Shear 

None 


Table IV. Environmental Properties for RD Test case: “Shelf Break”. 
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The vertical sound speed profile and bathymetry for this case are dis¬ 
played in Fig. 8(a). This characteristic is extended infinitely in the direction perpen¬ 
dicular to the page providing a 3-D shelf break environment. The locations of the 
source and two of the receivers are shown in Fig. 8(b), which is the same environment 
as Fig. 8(a) but viewed from the top. Isobaths are also drawn to provide a better 
perspective on the relative positions of the source and receivers. 



0 5 10 16 20 
Range (km) 


(a) 

Isobaths (m) 


400 350 300 250 200 150 100 100 



Range (km) 

(b) 

Figure 8. Sound speed profile: (a) Shelf break environment vertical view, (b) Shelf 
break environment horizontal view with sample locations of source and receivers. 
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d. Internal Waves 

This environment is defined by varying water column sound speed fea¬ 
tures and a flat, homogeneous bottom. More specifically this test case is the combi¬ 
nation of a simple, sinusoidally varying perturbation with a soliton wavefield, and is 
loosely based on the work of Tielburger, Finette, and Wolf [Ref. 32]. The aforemen¬ 
tioned environment is characterized by the following properties as shown in Table V. 


Bathymetry 

Flat bottom (zb = 200 m depth) 

Sound Speed Profiles 

varying 

Water Density 

1 g/cm 3 i 

Water Attenuation 

0.0 dB/A 

Bottom Sound-Speed 

1700 m/s 

Bottom Density 

1.5 g/cm 3 

Compressional Bottom 
Attenuation 

0.1 dB/km/Hz 

Shear 

None 


Table V. Environmental Properties for RD Test case: “Internal Waves”. 


The background (range-independent) sound speed profile used is the 
same as defined in Eq. (IV.35). The perturbations to this background are a com¬ 
bination of a sum of simple sinusoids and a train of soliton waves. The sinusoidal 
perturbations are defined by 

dcsin(z,r) = C (rg) e(~ 5 ) {cos (R'(i)r)} (IV.36) 

\BJ i=1 


where K(i) = [( 2 ooom)-(t^-i)( 300 r^ )I> B = 25m > and ^ is de fi ned suc h that the maximum 
value for dcsin is 7.5 m/s. 

The soliton perturbation is defined by 


dcsol(z,r) = C (4}) ^ 5Z |-^(*) 
\B/ *=i I 


seek 


( (R(i) -r) ' 
\ D(i) , 


(IV.37) 
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where B = 25 m, A(i) = lOe^ 0 ^" 1 )), D{%) = ^(“), J2(i) = fl(i - 1) - (7 - i)500 
for i = 2, • • •, 6, .R(l) = 14000m, and C is defined such that the maximum value for 
dc is 12.5m/s. 

The combination of the above perturbations is then simply defined by 

dc(z, r ) = dcsol(z, r) + dcsin(z, r ) (IV.38) 

and the result for the total sound speed structure is depicted in Fig. 9. Note that 
the solitons are modeled as propagating only in the plane of the page while the 
internal wave sinusoids are considered randomly oriented in all directions. For the 
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1620 
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1510 
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1500 

1495 

1490 


Figure 9. So un d speed profile: Combination of Sinusoidal and Soliton Perturbations. 


range-independent environments the source depth was between 140m and 160m and 
the locations of the receivers were at ranges between 2km and 8km and at depths 
between 40m and 60m. For the Shelf Break environment the source depth was at 
250m and the receivers were at a depth of 100m at ranges between 2km and 5km, 
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while for the Internal Waves environment the source depth was at 150m and the 
receivers were at a depth of 50m and at the same ranges as in the Shelf Break. Better 
insight to the scenarios is given in the next chapter, where the simulation results for 
the TDOA are presented. 
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V. TDOA SIMULATION RESULTS 


A. SYNTHETIC DATA 

In this section we present simulation results using synthetic data in the TDOA 
est ima tion problem. Our purpose is to distinguish how the two families of methods 
(classical - subspace) perform, what their basic differences are, and finally which 
method(s) produce the best overall results. Since the number of combinations using 
different signals, kind and amount of noise, and different characteristics of the signals 
and of the methods is huge, we have tried to provide some typical cases in order to 
provide insight into the behavior of the methods. We will try to explain how the 
different factors attribute to the final result, which is the estimation of TDOA. 

In the first subsection of the synthetic data we present cases from “Expo¬ 
nential”, “Sinusoidal”, “Damped Sinusoidal”, “Chirp” and finally “Damped Chirp” 
signals at different noise levels (SNR 15, 10, 5 dB). Our intent is to show the be¬ 
havior of the methods when the transient duration changes (short to long) and how 
successful they are when the desired TDOA to be predicted is large or small. In 
all cases there are different realizations of white noise with different variance added 
to each transient. For the subspace methods, the size of the covariance matrix (see 
Chapter 3) taken was 10 and 20 samples, and for the classical methods the number of 
segments (see Chapter 2) was either 1 or 2. In addition to the figures, corresponding 
tables s umm arize the results in order to provide a quick appreciation of each method’s 
behavior. 

In the second subsection of the synthetic data, we merge both kinds of methods 
(classical-subspace) to examine how the category of methods influences the other and, 
of course, if there is any substantial improvement in the results. 

1. Implementation Using Individual Methods 

The cases, whose results are presented in this subsection and in Appendix A, 
using synthetic data are the following: 
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• exponential transient with length of 100s and actual TDOA 50s; 

• exponential transient with length of 20s and actual TDOA 100s; 

• exponential transient with length of 50s and actual TDOA -10s; 

• sinusoidal transient with length of 100s and actual TDOA 50s; 

• damped sinusoidal transient with length of 100s and actual TDOA 50s; 

• chirp transient with length of 100s and actual TDOA 50s; and 

• damped chirp transient with length of 100s and actual TDOA 50s . 

Each case was tested with different realizations of additive noise and under three 
different levels of SNR. The figures and the comparative table for the first case is 
presented in this subsection. The remaining case results are submitted in Appendix 
A. Even though we have tested each kind of transient thoroughly changing various 
parameters, apart from the ones already defined, we have observed that the major 
influence to the desired result of TDOA is provoked by the noise (SNR), length 
of transient, amount of predicted TDOA, size of covariance matrix (for subspace 
methods), number of segments (for classical methods), and kind of transient. Having 
all the above in mind, we picked three typical cases from one kind of transient in order 
to see how the methods perform under different values of the aforementioned factors, 
and the same case from the rest of the transients for checking possible difficulties in 
estimating the TDOA for a particular type of signal. 

Figures 10 through 15 show TDOA results for the exponential transient for 
all methods using SNR values of 15dB, lOdB, and 5dB. The estimated values of time 
delay are fisted in Table VI. 
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Figure 10. Exponential transient: length L=100s; Actual TDOA=50s; SNR=15dB; 
White-Noise (cr%). (a) Subspace methods, Covariance Size=10. (b) Classical methods, 
number of segments=l. 
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Figure 11. Exponential transient: length L=100s; Actual TDOA=50s; SNR=15dB 
White-Noise (cr^). (a) Subspace methods, Covariance Size=20. (b) Classical methods 
number of segments=2. 













Figure 12. Exponential transient: length L=l00s; Actual TDOA=50s; SNR=10dB; 
White-Noise (cr%). (a) Subspace methods, Covariance Size=10 (b) Classical methods, 
number of segments=l 
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Figure 13. Exponential transient: length L=100s; Actual TDOA=50s; SNR=10dB; 
White-Noise (<j l). (a) Subspace methods, Covariance Size=20. (b) Classical methods, 
number of segments=2. 














































Figure 14. Exponential transient: length L—100s; Actual TDOA=50s; SNR=05dB; 
White-Noise (cr%). (a) Subspace methods, Covariance Size=10 (b) Classical methods, 
number of segments=l 
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Figure 15. Exponential transient: length L=100s; Actual TDOA=50s; SNR=05dB; 
White-Noise (cr^). (a) Subspace methods, Covariance Size=20 (b) Classical methods, 
number of segments=2 














































Method 

Parameters 

SNR 15 dB 

SNR 10 dB 

SNR 5 dB 

MUSIC 

Cov-Mat: 10 

50 

51 

51 

Cov-Mat: 20 

50" 

49 

49 

Modified MUSIC 

Cov-Mat: 10 

50 

51 

52 

Cov-Mat: 20 

50 

49 

50 

MIN-NORM 

Cov-Mat: 10 

50 

51 

50 

Cov-Mat: 20 

50 

49 

49 

PCLP 

Cov-Mat: 10 

50 

51 

52 

Cov-Mat: 20 

50 

49 

49 

ESPRIT 

Cov-Mat: 10 

49 

51 

50 

Cov-Mat: 20 

50 

50 

50 

Root-MUSIC 

Cov-Mat: 10 

50 

50 

49 

Cov-Mat: 20 

50 

49 

49 

SCOT 

#-Segs: 1 

50 

50 

92 

#-Segs: 2 

50 

52 

47 

PHAT 

#-Segs: 1 

50 

50 

92 

#-Segs: 2 

50 

52 

47 

Roth 

#-Segs: 1 

-299 

-502 

248 

#-Segs: 2 

50 

52 

-204 

X-Corr 

#-Segs: 1 

50 

50 

53 

#-Segs: 2 

50 

50 

47 


Table VI. Exponential transient: length L=100s; white noise (cr^ = 2); actual 
TDOA=50s. 

After examining the corresponding Figs. 10 - 15 and Table VI for the first case 
of the synthetic data we can say that almost all the methods performed satisfactorily. 
The figures for the subspace methods in all cases, Figs. 10(a) - 15(a), were clear 
without any subsidiary peaks, compared with the respective figures of the classical 
methods, Figs. 10(b) - 15(b). For medium noise level (SNR=10dB), Figs. 12 - 13, 
and even more for high noise level (SNR=5dB), Figs. 14 - 15, the subspace methods 
gave more accurate results. The increase of the size of the covariance matrix had a 
slight improvement in the performance for the subspace methods. On the contrary 
the increase of the number of segments for the classical methods had more positive 
influence on their performance, especially when the SNR became very low (5dB). Of 
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the classical methods, Roth had the worst performance and cross-correlation had the 
best performance, on average. Subspace methods with both approaches had similar 
behavior in all conditions. 

In summary, after reviewing all cases of the synthetic data (see Figs. 51 - 
86 and Tables XX - XXIV in Appendix A) for the performance of all methods, the 
following observations can be made. 

• The subspace methods uniformly had very efficient performance compared 
with the classical methods 

• The subspace methods were more consistent in the quality of the results as 
a group. The classical methods, on the other hand, in some cases had great 
diversity among them in the quality of their performance. Among the clas¬ 
sical methods, Roth had the poorest performance; and in general the cross 
correlation was more consistent in accuracy than the other methods. 

• The plots of the subspace methods were more distinct and it was easier to iden¬ 
tify the peak corresponding to the correct Time Difference of Arrival (TDOA). 

• In cases where both methods provided the expected results the classical meth¬ 
ods had more accurate peaks than the subspace methods. However both meth¬ 
ods provided results in which accuracy was completely acceptable. 

• In order to perform with multiple segments the classical methods needed to 
have appropriate length of input data (signal plus noise). Otherwise phe¬ 
nomena of aliasing appeared especially in cases with a large “value” of TDOA 
between the two data sets. This drawback did not affect the subspace methods 
at all since they do not deal with segmentation. 

• Methods like root MUSIC or ESPRIT can produce direct numerical estimates 
for the time delay without the need to search for a peak. 

• Subspace methods were more “resistant” to noise compared to the classical 
methods. 

• For short duration of transients the performance of the methods decreased in 
accuracy, especially when it was combined with low SNR. In this case classical 
methods gave totally wrong results compared with the subspace methods, 
whose performance was inside the permissible limi ts 0 f error (Table XX). 

• For the case of a small value of predicted TDOA (Table XXI) the subspace 
methods performed more than satisfactorily even under high “noise condi¬ 
tions”. The performance of the classical methods was similar with that of 
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the subspace methods for high SNR, but they became quite inferior when the 
noise was increased. 

• In general the increase in the size of the covariance matrix in subspace methods 
improved their performance, especially when the SNR was small. Better results 
were achieved when the size of the covariance matrix was greater than 20, but 
increased the total number of computations something that made the subspace 
methods to delay enough in order to provide us with the desired result. 

For the cases where all the remaining factors are the same (except for the type of 
transient), Tables XXII, XXIII, XXIV, and Figs. 10 - 15 and 63 - 86, we can say 
that all the methods have almost the same behavior with different transients as with 
the exponential transient. Only the sinusoidal transient had noticeably worse perfor¬ 
mance than the others and was influenced more with high noise level than the rest of 
the transients. 

2. Implementation Using Combination of Methods 

Before we continue with the MMPE data, it is worth examining if the “com¬ 
bination” of the two fa mili es of methods, subspace and classical, perform better than 
each one separately. In our case the word “combination” means that before we apply 
the subspace methods in the frequency domain (cross-spectrum), the desired weight¬ 
ing factor will be implemented according to the processor we would like to use. The 
main idea was to examine if the subspace methods, with the assistance of the weight¬ 
ing factor, could provide us with more accurate results (peaks). Several test cases 
were conducted and it was determined that the peaks were more accurate but the 
results were not better them without implementing the combination. As a general 
remark, each weighting factor influenced the performance of the subspace methods 
differently in such a way that the subspace methods lost their good consistency in 
results. Figures 16 - 17 and the corresponding Table VII give a sample of the “com¬ 
bination” performance. 
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Method 

SCOT 

PHAT 

ROTH 

CROSS 

CORRELATION 

NONE 

MUSIC 

95 

95 

-879 

103 

103 

Modified-MUSIC 

96 

96 

-600 

103 

103 

MIN-NORM 

90 

90 

-620 

103 

103 | 

PCLP 

91 

91 

-871 

103 

103 

ESPRIT 

95 

95 

-728 

103 

103 

Root-MUSIC 

96 

96 

-864 

103 

103 


Table VII. Damped Sinusoidal transient: length L=400s; white noise (cr 2 = 2); 
SNR=10dB; actual TDOA=100s. 
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Figure 16. Damped Sinusoidal transient: Actual TDOA=100s; SNR=10dB; White- 
Noise (cr„) Covariance Size=10, number of segments=l. (a) Subspace methods with 
SCOT, (b) Subspace methods with PHAT. 
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Figure 17. Damped Sinusoidal transient: Actual TDOA=100s; SNR=10dB; White- 
Noise (al) Covariance Size=10, number of segments=l. (a) Subspace methods with 
ROTH, (b) Subspace methods with Cross Correlation. 
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B. MMPE DATA 

In this section results from data generated by the MMPE propagation model 
are presented. The particular environments that were used to generate these data 
are the ones presented in Chapter 4. These simulations were conducted to verify 
the observations from the synthetic data and to test the methods with broadband 
transients under more complicated environmental conditions. For each environment, 
we provide a Transmission Loss (TL) plot relative to the depth and the frequency 
bandwidth, the waveforms of the source and of the two receivers in the time domain, 
and the results of all TDOA methods. It will be seen that in all cases the signals 
in each location have totally different shape and magnitude due to the multipath 
structure of the waveguide. As mentioned before, the environments are divided into 
two categories: range-independent (the received signal depends only on the distance 
between the source and the receiver [Flat - Sound Channel cases], and range-dependent 
(the sound speed profile changes according to range and direction from the source 
[Shelf Break - Internal Waves cases]. Two cases are considered for each environment. 
For the range-independent environments the receivers axe chosen to be at different 
ranges from the source. For the range-dependent environments the receivers are 
selected to be at the same ranges but in different relative directions from the source. 
In this way the difference between the two categories will be more evident. For more 
detailed description of the MMPE environments see Chapter 4. 

1. Flat Bottom 

Two cases are considered for this environment. The receivers are at distances 
of 2km and 3km, and 2km and 5km for each case, respectively. In all cases different 
realizations of white noise were added to each transient with SNR=10dB. In Figs. 18 
and 19 are shown the Transmission Loss (TL) plot of depth versus frequency and the 
waveforms at the source and at the two receivers in the time domain for one of the 
two cases. The results for both cases are shown in Figs. 20 and 21 and summarized 
in Table VIII. From these results, it appears that all methods perform in a similar 
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manner to that in the section with the synthetic data. All of them gave results very 
close to the actual value for both cases with slightly better accuracy by the subspace 
methods. The difference in the figures between the two group of methods is still 
the same. Subspace methods have more clear plots and it is easier to pick up the 
correct peak from their graphs than the corresponding ones from classical methods 
(subsidiary peaks.) In this environment the root-finding subspace methods (ESPRIT 
- rootMUSIC) had the best performance of all the methods and Roth gave the worst 
results compared with the rest. 


Method 

Estimated TDOA 

Estimated TDOA 

Case 1 

Case 2 

ACTUAL 

0.6936 

2.1209 

MUSIC 

0.6953 

2.1250 

Modified-MUSIC 

0.6641 

2.0879 || 

MIN-NORM 

0.6943 

2.1240 

PCLP 

0.6943 

2.1240 

ESPRIT 

0.6934 

2.1221 

Root-MUSIC 

0.6943 

2.1230 

SCOT 

0.7324 

2.1182 

PHAT 

0.7324 

2.1182 

Roth 

0.5713 

2.1328 

X-Corr 

0.7285 

2.1133 1 


Table VIII. TDOA results for Flat Bottom isospeed case. 
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Figure 18. Transmission Loss for FLAT Bottom environment. 
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Figure 19. Time domain waveforms for FLAT Bottom environment, (a) Source, (b) 
Receiver 1. (c) Receiver 2. 
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Figure 21. FLAT Bottom environment Subcase 2: Actu al TDOA=2.1209s; 
SNR=10dB; White-Noise (<r%). (a) Subspace methods, Covariance Size=10. (b) Clas¬ 
sical methods, number of segments=l. 













2. Sound Channel 

Two cases are considered for this environment. In the first case the receivers 
are at distances 2km and 4km, while in the second case they are at 2km and 8km. In 
all cases a different realization of white noise was added to each transient to obtain 
a SNR of lOdB. Figures 22 and 23 show the Transmission Loss (TL) plot of depth 
versus frequency and the waveforms of the source and of the two receivers in the 
time domain for one of the two cases. The results for both cases are shown in Figs. 
24 and 25 and s umm arized in Table IX. The results show us a similar picture for 
the performance of the methods. The only difference is that in this environment the 
“gap” between the estimated values of TDOA from the classical methods and the 
actual predicted value becomes larger. 


Method 

Estimated TDOA 

Estimated TDOA 

Case 1 

Case 2 

ACTUAL 

1.5830 

4.5929 

MUSIC 

1.5830 

4.5967 

ModifiedL-MUSIC 

1.5781 

4.5967 

MIN-NORM 

1.5830 

4.5967 

PCLP 

1.5830 

4.5967 

ESPRIT 

1.5820 

4.5967 

Root-MUSIC 

1.5820 

4.5967 

SCOT 

1.5771 

4.5527 

PHAT 

1.5571 

4.5527 

Roth 

1.5596 

4.5527 

X-Corr 

1.5732 

4.5527 


Table IX. TDOA results for SOUND CHANNEL case. 
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Figure 22. Transmission Loss for SOUND CHANNEL environment. 
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Figure 23. Time domain waveforms for SOUND CHANNEL environment, (a) Source, 
(b) Receiver 1. (c) Receiver 2. 
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Music: 1.583 


MinNorm: 1.583 


M Music: 1.5781 


-3 *2 -1 



-3 -2 -1 


Figure 24. SOUND CHANNEL environment Subcase 1: Actual TDOA=1.5830s; 
SNR=10dB; White-Noise (cr*). (a) Subspace methods, Covariance Size=10. (b) 
Classical methods, number of segments=l. 
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Figure 25. SOUND CHANNEL environment Subcase 2: Actual TDOA=4.5967s; 
SNR=10dB; White-Noise (al). (a) Subspace methods, Covariance Size=10. (b) 

Classical methods, number of segments=l. 





























3. Shelf Break 

For this environment the receivers are at distances 2km and 3km but in posi¬ 
tions relative to the source as shown in Fig. 26 for each of the two cases considered. 
In particular Fig. 26 provides us with a top view of the geometry for the two cases, 
indicating the positions of the source and of the receivers relative to the source, noting 
also the corresponding depths and the respective ranges for each case. Once again, 
different realizations of white noise were added to each transient to obtain a SNR 
of lOdB. Figures 27 and 28 show the Transmission Loss (TL) plot of depth versus 
frequency and the waveforms of the source and of the two receivers in the time do¬ 
main for one of the two cases. The results for both cases are shown in Figs. 29 and 
30 and s umm arized in Table X. As we now get into more complex environments 
(range-dependent), it is obvious that the graphs of the classical methods become 
more “noisy” than the ones in the range-independent enviro nm ents, and also that 
they had an abrupt decrease in their performance except for the Cross-correlation 
method. Subspace methods still gave us accurate results and all of them were very 
consistent about the small diversity of their estimated values of TDOA. 


LEGEND _ 

Case 1 Receivers (§) 
Source <£> 

Case 2 Receivers<£) 
Source Depth : 250m 

Receiver Depth: 100m 



Figure 26. Top view of positions of source and receivers for both cases of S HE LF 
BREAK environment. 
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Depth (m) 


Source Level (cB re 1 m) 



340 360 393 400 420 440 460 


Freq (Hz) 

Figure 27. Transmission Loss for SHELF BREAK environment. 


Method 

Estimated TDOA 

Estimated TDOA 

Case 1 

Case 2 

ACTUAL 

0.7173 

0.7407 

MUSIC 

0.7139 

0.7402 

Modified-MUSIC 

0.7139 

0.7490 

MIN-NORM 

0.7139 


PCLP 

0.7139 


ESPRIT 

0.7139 

0.7471 

Root-MUSIC 

0.7139 

0.7471 

SCOT 


-0.2793 

PHAT 


-0.2793 

Roth 

1.2998 

1.6211 

X-Corr 

0.8154 

0.7803 


Table X. TDOA results for SHELF BREAK case. 
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(c) 

Figure 28. Time domain waveforms for SHELF BREAK environment, (a) Source, 
(b) Receiver 1. (c) Receiver 2. 
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sical methods, number of segments=l. 
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4. Internal Waves 

For this environment both receivers are at a distance of 3km but at different 
positions relative to the source as shown in Fig 31 for each of the two cases considered. 
Figure 31 provides a top view of the geometry for the two cases, indicating the 
positions of the source and of the receivers relative to the source, noting also the 
corresponding depths and the respective ranges for each case. As with the previous 
environments, different realizations of white noise were added to each transient to 
obtain a SNR of lOdB. Figures 32 and 33 show the Transmission Loss (TL) plot of 
depth versus frequency and the waveforms at the source and at the two receivers in 
the time domain for one of the two cases. The results for both cases are shown in Figs. 
34 and 35 and s umm arized in Table XI. Our observations for this environment are 
almost the same as with the Shelf Break. The only difference is that even the subspace 
methods give slightly worse results compared with the actual value of TDOA. This 
was expected since the Internal Waves was the most complex environment used in 
this thesis. These results are still within the permissible limits of error. 



Case 1 Receivers (§) 

Source 

O 

Case 2 Receivers (£) 

Source Depth : 

150m 

Receiver Depth: 

50m 


Figure 31. Top view of positions of source and receivers for both cases of INTERNAL 
WAVES environment. 
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Depth (m) 


Source Level (cB re 1 m) 



340 360 380 400 420 440 460 

Freq (Hz) 


Figure 32. Transmission Loss for INTERNAL WAVES environment. 


Method 

Estimated TDOA 

Estimated TDOA 

Case 1 

Case 2 

ACTUAL 

0.7770 

0.2109 

MUSIC 

0.7734 

0.2090 

Modified-MUSIC 

0.8311 

0.1885 

MIN-NORM 

0.7822 

0.2061 

PCLP 

0.7793 


ESPRIT 

0.7959 


Root-MUSIC 

0.7930 

0.2041 

SCOT 

hheeshhi 

0.1289 

PHAT 


0.1289 

Roth 

2.5938 


X-Corr 

0.8574 

0.1357 


Table XI. TDOA results for INTERNAL WAVES case. 
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(C) 

Figure 33. Time domain waveforms for INTERNAL WAVES environment, (a) 
Source, (b) Receiver 1. (c) Receiver 2. 
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0 


MinNorm: 0.78223 


M Music: 0,83105 



Figure 34. INTERNAL WAVES environment Subcase 1: Actual TDOA=0.7770s; 
SNR=10dB; White-Noise (cr%). (a) Subspace methods, Covariance Size=10. (b) 
Classical methods, number of segments=l. 





























































































































Figure 35. INTERNAL WAVES environment Subcase 2: Actual TDOA=0.2109s; 
SNR=10dB; White-Noise (a%). (a) Subspace methods, Covariance Size=10. (b) 
Classical methods, number of segments=l. 










































































Taking into consideration the figures and the tables with the TDOA results 
from all the environments, we observe that all the methods performed quite accurately. 
It should be noted that the subspace methods were more consistent than the classical 
ones. This can be verified from the range-dependent environments, where all the 
subspace methods gave results more close to the actual value and they had only a small 
dispersion among them. On the contrary, the classical methods had slightly inferior 
performance in the range-independent environments and much worse in the range- 
dependent ones, with a wider dispersion of values. This difference in the performance 
among all the methods between the two categories of the environments also appears 
in the results of the localization problem in the next chapter. 
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VI. THE LOCALIZATION PROBLEM 


In this chapter, we use the TDOAs between the receivers, estimated by the 
various methods, in order to estimate the location of the target. For this purpose a 
homonymous algorithm will be implemented whose main characteristic is to provide 
us the passive assessment of target position by using a corresponding set of equations, 
called Time of Arrival (TOA) equations. The TOA equations try to correlate the 
positions of the receivers with the position of the target using the time that a transient 
needs to travel from one location to the other. In other words, these equations try 
to fin d a correspondence between the signal arrival times and the target range. By 
solving the TOA equations for only two receivers, we are not going to get one specific 
location of the target, but an infinite number of possible locations since this solution 
corresponds to a bearing from this particular pair of receivers. Under these conditions, 
in order to e limin ate this ambiguity we need to solve the TOA equations for all the 
receivers at the same time, so the desired target position will be the result of the 
crossing of all the bearings. This expectation is quite natural since, as far as we 
know, only one solution (target position) can account simultaneously for the travel 
time of the transient to each receiver. The algorithm that is implemented is simpler 
than most since it transforms the original nonlinear TOA equations into a set of linear 
equations. The inputs (receiver position and signal’s arrival time for each receiver) 
correspond to a specific set of outputs (target position and transmission time of the 
tr ansi ent) - The rest of this chapter is based on a former thesis [Ref. 33] related to 
this subject. 

A. PRESENTATION OF TDOA ALGORITHM 

In this section, a brief description of the TDOA algorithm will be provided to 
the reader since a detailed derivation of this algorithm is not within the scope of this 
thesis. Before we continue with the rest of our discussion, it would be advisable to 
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make clear that the solution to the localization problem, which is attempted to be 
given here, is conducted in two dimensions in the x - y plane. Let us assume that 
we have one source (target) and N receivers (buoys). The location for each receiver 
is given by the vector 

(VI-1) 

where x { and y* are the x-y coordinates of buoy i. Let the arrival time of the transient 
at buoy i, defined as the time it takes the leading edge of the transient to travel from 
the source to the receiver, be denoted as In a similar manner, the transmission 
time of the transient from the target will be denoted as t, and the desired position of 
the target is defined by 

(VI-2) 

where x equals the target’s x coordinate position, y equals the target’s y coor dina te 
position. 

In order to estimate the desired variables, t and r, we have to solve the following 
equation in matrix form 

Ar = qf + s, (VI.3) 

where 

-x 2 yi - y 2 
x 2 2/2 — 2/2 

x n -i — x n 2/jv-i — lI n 
ti-t 2 

h ~ h 

tN-l ~ tiV 


(VI.5) 



(VI-4) 
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and 



INI 2 -INI 2 -A5 + A§ 
INP-IhlP-Ai + c 2 ^ 


(VI.6) 


Lii^ip-imp-a^+^j 

The variable c represents the sound speed, which in our case will be assumed constant 
throughout the whole computational grid. Equation (VL3) is simply the transforma¬ 
tion of the nonlinear TOA equation 


(Ii - X? + ( Vi - yf - c 2 (t, - t) 2 = 0 (VI.7) 


into a linear TDOA equation which relates target position to transmission time 

x{x n X n +i )~\~y(yn Vn+ 1) = C ^n+l) "b fan '^•'71+1 "b Vn Vn+1 ^ ^n+ 1)* 

(VI.8) 

T his equation is formed by subtracting two equations like (VI.7) for successive pairs 
of the N receivers. 

The solution of Eq. (VL3) is computed in a least squares sense in order to 
find the target range r in terms of the transmission time t. Equation (VI.3) becomes 

r = g£ + h, (VI.9) 

where 

g = (A) + q. (VI.10) 

h = (A)+s, (VI.11) 

and + denotes the pseudoinverse. Here we make the distinction that when N = 3 
the expected solution is exact for the corresponding set of equations, but for the case 
N > 3 the result is the least squares solution, which minimizes the equation error for 
the same set of Eqs. (VI.7). 

After some calculations, we come up with the following formula for the k th 
range equation: 

ll r k — gt — h| | 3 = c s (tk — t) s , (VI-12) 
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which after expanding and simplification produces 

at 2 + 2bt + d = 0 , 

(VI. 13) 

where 

<* = ||g|| 2 - c 2 , 

(VI.14) 

b = g T • h - g T • r k + c 2 t k , 

(VI. 15) 

d— ||h — r k || 2 - c 2 t k . 

(VI. 16) 


By solving the quadratic formula (VI. 13), we find two possible solutions for transmis¬ 
sion time t, but only one of them is correct. We must use the correct solution in Eq. 
(VL9) in order to find the location of the target. 

In theory, the solution for t is independent of which particular range equation 
(k) was used to develop the quadratic equation. In practice, because of measurement 
errors, the solutions will be different. Therefore, it is advisable to develop N - 1 
estimates for t using k = 1,2, • • •, N — 1, discarding any erroneous values and average 
the remaining solutions. This multiple solution technique also resolves which root is 
correct. The result can then be used in (VL9) to find target position. 

B. MATLAB IMPLEMENTATION 

This section explains how the localization testing was implemented in MAT- 
LAB using data produced by the MMPE propagation model. The problem was rep¬ 
resented by a “target”, which is the sound source, and two or more sonobuoys, which 
are the receivers. More specifically, we specified a cartesian coordinate system and 
set the true position of the target. The buoy positions were given to the MATLAB 
program relative to the target’s position. This was done by having as inputs the 
quadrant number, the distance (range) to the buoy, and the distance in x-coordinate 
to the buoy. In this way, it was possible to create any geometric scenario in two 
dimensions. 
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After specifying the geometry of the scenario, we had to find the value of 
the sound speed that would be used for each environment. To be as realistic as 
possible with each case, we evaluated the group speed of a specific mode that would 
dominate during propagation of the signal. This was executed for each water column 
corresponding to each buoy location in the range-dependent cases and only once in 
range-independent cases, since the sound speed profile remains the same throughout 
the computational grid. For the range-dependent environments, we averaged those 
values of the group speed for each buoy in order to end up with one value to use for 
the rest of the algorithm. 

After these steps, the data from MMPE (after being modulated/demodulated 
as described in Chapter 4) was read into the localization program. Then using all 
the previously described methods (subspace - classical) we evaluated the TDOAs for 
all combinations of pairs of buoys. Since we don’t know absolute TOA for any of the 
buoys, we set TOA for the 1st buoy equal to zero, TOA( 1) = 0, and measured the 
relative TOA for the other buoys with respect to the first buoy as a reference. 

The rest of the algorithm continues as it was presented in the last section to 
evaluate the “Estimated Target Position.” The entire procedure is repeated for all the 
TDOA methods for each environment and for each case and using the same signals. 
In this way, it is possible to evaluate the performance of all the methods under the 
same circumstances and compare results among them. Figure 36 summarizes all the 
steps of the localization procedure. 

C. SIMULATION RESULTS 

In this section, simulation results for the localization problem using the four 
MMPE environments are presented. For each environment, there are two cases con¬ 
sidered. In the first, the minimum number of three sonobuoys is used. (Recall that 
three is the minimum , since we require at least two separate TDOA measurements to 
locate the target.) In the second case, more than three sonobuoys are used. 
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Create 


Find Average 


Evaluate TDOA 

Geometric 

- p 

Group Speed 


-- p 

For all pair 

Scenario 


for Buov Locations 



combinations 

1 



Figure 36. Flowchart of Localization Algorithm implemented in MATLAB environ¬ 
ment using data from MMPE. 

For each environment and each TDOA method we show a pair of figures with 
the locations of the buoys, the real target location and the estimated location from 
each method. Each figure in the pair corresponds to the two cases mentioned above (3 
buoys and more than 3 buoys), and contains two subfigures with the results from the 
subspace and from the classical methods respectively. Each environment is accompa¬ 
nied by a table. The tables contain the error in target’s position, i.e., the diff erence 
between the actual location and the estimated one for each method. In this way there 
is a fair evaluation for all the methods, since for this kind of problem the difference 
of the two positions (real-estimated) is of most concern. 
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1. Flat Bottom 

The characteristics of this environment are described in Chapter 4. The true 
geometric positions for two scenarios with 3 and 5 receivers (buoys) are presented in 
Figs. 37 and 38 respectively. The source was at a depth of 150m while the receivers 
were at a depth of 50m and ranges of 2, 3 and 5 km from the source. Table XII has 
the target position error for all methods for both cases in this environment. 

Observing Table XII, we conclude that all methods performed efficiently, but 
the position accuracy was better for the subspace methods than the classical meth¬ 
ods. In both groups of methods there was an improvement with the increase of the 
number of sonobuoys, except from the Roth method. The variance of the results from 
the subspace methods was much smaller than the respective one from the classical 
methods. Further, the cross-correlation (X-Corr) had the best performance from the 
classical methods and was the only one from this family whose results were closer to 
those of the subspace methods. 


Method 

Target Position Error in (m) 


3 buoys 

5 buoys 

MUSIC 

29.68 

23.74 

Modified-MUSIC 

47.46 

24.67 

MIN-NORM 

30.16 

22.45 

PCLP 

30.16 

22.42 

ESPRIT 

30.76 


Root-MUSIC 

30.06 

21.81 

SCOT 

100.87 

67.52 

PHAT 

100.87 

67.52 

Roth 

87.17 

183.81 

X-Corr 

37.04 

31.42 


Table XII. Target position Error for Flat bottom isospeed case. 















































(w) A 


1st buoy 
2nd buoy 
3rd buoy 
4th buoy 
5th buoy 
Target posit ton 
Est-Pos-MUSIC 
Est-PosMod MUSIC 

Est-Pos MN 
Est-Pos-PCLP 
Est-Pos-ESP 
Est-Pos-rootMU 


-10000 l— 
-10000 



— 

— 

— 

— 

— 

O 1st buoy 
+ 2nd buoy 
x 3rd buoy 







A 4th buoy 

0 6th buoy 

O Target posit ton 







-f Est-Pos-SCOT 
+ EstPosPHAT 
a Est-Pos-ROTH 
ts EstPos-XCORR . 



$ 








— 

.O ■ 







-£.- 

. 

_ 

_*_ 






_ 






1- 

i - i 

I_ 

i_ 

i_ 

i_ 

1_ 


Figure 38. Localization problem for Flat bottom isospeed case with 5 sonobuoys. (a) 
Subspace methods, (b) Classical methods. 
























2. Sound Channel 

The characteristics of this environment are described again in Chapter 4. The 
true geometric positions of the two scenarios with 3 and 5 receivers (buoys) are shown 
in Figs. 39 and 40, respectively. The source was at a depth of 150m and the receivers 
were at a depth of 50m and at ranges of 2 and 4 km from the source. The results are 
also accompanied by Table XIII presenting the Target Position Error for all methods. 
Once again the subspace methods had better accuracy than the classical methods 
as a group. In both families of methods, we observe that increasing the number of 
the buoys does not, by itself, assist the algorithm to more accurately predict the 
target position. The main consideration that has to be made is the disposition of 
the receivers and then the number of them. If the disposition is suitable, then an 
increase of the number of the sonobuoys may improve the result since, as we have 
seen in section A of this chapter, Eq. VI.3 will be solved either uniquely if N = 3 or 
in a least squares sense if N > 3. The performance was similar to the Flat bottom 
case, with Roth being more unstable in the results, and SCOT, PHAT performing 
better than X-Corr. 


Method 

Target Position Error in (m) ]j 


3 buoys 

5 buoys 

MUSIC 

19.65 

25.72 

Modified-MUSIC 

21.35 

28.42 ! 

MIN-NORM 

19.65 

26.65 ! 

PCLP 

19.65 

26.65 

ESPRIT 

20.50 

26.20 

Root-MUSIC 

19.65 

26.65 

SCOT 

47.04 

53.02 

PHAT 

47.04 

53.02 

Roth 

24.70 

487.13 

X-Corr 

45.92 

65.09 


Table XIII. Target position Error for Sound Channel case . 
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Figure 39. Localization problem for Sound Channel case with 3 sonobuoys. (a) 
Subspace methods, (b) Classical methods. 


— 

— 

— 

— 

— 

_ 

— 

o 1st buoy 
+ 2nd buoy 
x 3rd buoy 

O Target posit ton 
+ Est-Pos-SCOT 
+ Est-Pos-PHAT 

A Est-Pos-ROTH 

0 Est-Pos-XCORR 







— 




O 








♦ 

— 







X 











_ 





_ 

_ 





i_ 

i_ 

1_ 

i_ 



1_ 


91 























































3. Shelf Break 

Using the characteristics for this environment discussed in Chapter 4, the true 
geometric positions of two scenarios with 3 and 6 receivers (buoys) are presented in 
Figs. 41 and 42 respectively. The source was at a depth of 250m and the receivers 
were at a depth of 50m and at range of 3 km from the source, located upslope and 
downslope concurrently (see similar Fig. 26). Table XIV summarizes the results of 
the performance for all methods. In this range-dependent environment, we see that 
subspace methods have similar accuracy in their results compared with the one in 
the range-independent environments. In comparison with the subspace methods, the 
classical methods performed much poorer. This could be expected if we recall their 
quality of results in the estimation of TDOA (see Table X.) In general, there was an 
improvement of the results for all methods (except X-Corr) with an increase in the 
number of the sonobuoys. This improvement was more drastic for Roth, something 
that also was expected considering the diversity of the results from this method. An¬ 
other remark is that because of the complexity of the range-dependent environments 
the buoys disposition is very critical. An inappropriate disposition can create larger 
errors in target position than the range-independent environments. 


Method 

Target Position Error in (m) 


3 buoys 

6 buoys 

MUSIC 

24.65 

22.48 

Modified-MUSIC 

31.04 

22.64 

MIN-NORM 

25.32 

20.45 

PCLP 

25.32 

20.48 

ESPRIT 

27.29 

21.96 

Root-MUSIC 

27.29 

21.89 

SCOT 

434.35 

432.77 

PHAT 

434.35 

432.77 

Roth 

787.43 

418.26 | 

X-Corr 

93.53 

215.67 


Table XIV. Target position Error for Shelf Break case . 
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Figure 41. Localization problem for Shelf Break case with 3 sonobuoys. (a) Subspace 
methods, (b) Classical methods. 
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4. Internal Waves 

Using the characteristics for this environment described in Chapter 4, the true 
geometric positions of two scenarios are shown with 3 and 6 receivers (buoys) in 
Figs. 43 and 44, respectively. The source was at a depth of 150m and the receivers 
were at a depth of 50m and at range of 2 km from the source located on both the 
sides. We have to mention that some of the receivers were located in the section 
of the computational grid with sinusoidal and soliton perturbations and the rest 
of them were in the other section with sinusoidal perturbations only (see similar 
Fig. 31). Table XV contains the results of the performance for all methods. This 
environment is the most complex that we used so far and it deviates considerably from 
the basic concept of the localization algorithm, which assumes a constant value for the 
sound speed. For this environment all the methods have decreased performance. The 
difference between subspace and classical was the same as before, shown in Table XI. 
The only difference is that the increase in the number of the sonobuoys with different 
disposition had worse results than in the Shelf Break environment and that X-Corr 
gave s imil ar results to the subspace methods. 


Method 

Target Position Error in (m) || 


3 buoys 

6 buoys 

MUSIC 

65.90 

160.52 

Modified-MUSIC 

64.64 

153.58 

MIN-NORM 

63.02 

164.64 

PCLP 

62.29 

165.85 

ESPRIT 

63.81 

164.78 

Root-MUSIC 

66.11 

■E1E3HII 

SCOT 

423.57 

849.08 

PHAT 

423.57 

849.08 

Roth 

226.81 

515.59 

X-Corr 

49.22 

160.37 |i 


Table XV. Target position Error for Internal Waves case . 
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Figure 44. Localization problem for Internal Waves case with 6 sonobuoys. (a) 
Subspace methods, (b) Classical methods. 




























Prom the corresponding figures and tables for each environment, one can con¬ 
clude that in most cases all the methods performed quite successfully. The subspace 
methods were more consistent in their results as a group; this is something that didn’t 
happen with the classical methods. Further, with different realizations of the noise, 
the variance of the results of the subspace methods was much smaller than for the 
classical methods. The location accuracy in the range-dependent environments was 
less than in the range-independent environments and was reduced even more when 
the distance between receivers and source became greater (4 to 8km). Another major 
factor for the successful localization of the target was the geometrical disposition of 
the receivers relative to the source position. In other words, if the receivers were 
located in such a way that the bearing from each pair formed an angle of between 
60° and 120°, then the estimate of their crossing will have less error than otherwise. 

A larger number of receivers does not always contribute to more accurate 
detection. T his can be justified from the fact that when N is equal to 3, the solution 
from Eq. (VI.9) is unique, while a value of N > 3 requires a least squares solution 
to minimize the equation error in the over specified set of equations. As discussed in 
Section A of this chapter, the solution for t is independent of which particular range 
equation (k) is used to develop the quadratic equation, and in order to overcome 
the obstacle of measurement errors, we have to develop N — 1 estimates for t using 
k = 1,2,**-, TV" — 1, discarding any erroneous values and averaging the remaining 
solutions. A larger value of N does not always lead to a better result. In other words 
for more accurate results, the placement of the buoys is more important than the 
number of the buoys. 
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VII. TRACKING PROBLEM 


A. DOPPLER IMPLEMENTATION IN MMPE 

In actual underwater acoustics problems, the source and/or the receiver are 
not stationary, but they are in motion. This motion causes temporal fluctuations, 
something we have to take into consideration if we want to find out how the under¬ 
water acoustic channel influences various applications such as communications [Ref. 
34]. In our case we will consider motion only due to the source (target), and the 
receivers (buoys) will be stationary. Figure 45 gives a characteristic example of the 
effect. The result of the source motion is to provoke a change in the value of the actual 
source frequency observed, /$(#), in other words in the received frequency at buoy i 
compared with the original transmitted frequency, fy. The received frequency, fs(Q), 
will depend on the transmitted frequency, fa, the observation angle, 6, the source 
speed, vs, and the direction of motion, <f>s [Ref- 35], 


fs{0) = 


_ fr _ 

1 - (Jj) cos(6 - <ps) 


(vn.i) 


I 



Figure 45. Effect of Doppler due to linear motion of the source. 
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According to [Ref. 35] when |0 S | < 90° then the observed frequency bandwidth 
downrange, BW T S + , will be equal to 

BW S + = BWt + ^ (f T ,max + f T ,min |sin0 S |), (VII.2) 

where BW T is the actual transmitted bandwidth, and f T ,max and f TtJnin are the actual 
transmitted maximum and minimum frequencies, respectively, of the tr ans mitted 
band. The observed center frequency downrange is defined as [Ref. 35] 

fs,c = fr,c + y c (/r, max fTymin I sin 0s|), (VII.3) 

where f T , c is the actual transmitted center frequency. When |0 5 | > 90°, then the 
observed bandwidth downrange will be [Ref. 35] 

BW r s ~ = BWt + Y (f T ,max I sin 4) S | 4- , (VII.4) 

and the observed center frequency is equal to 

Us 

m = fr,c + — {fr,max I Sin 0s | + f T ,min) ■ (VII. 5) 

The equation for the PE starting field in the vertical wavenumber domain has 
the general form [Ref. 35] 


1> ( r = 0, /) = e~ ikzZ3/ ip 0 (k z , /) - e ikzZ, ijj 0 (k z , /), (VII.6) 


where the functions ip 0 (k z , f ) represent the starting field in free space. The role of the 
exponentials in Eq. (VII.6) is the creation of the appropriate interference structure 
between the source (target), which lies at depth z s and its image about the plane of 
the free surface (z — 0). The vertical wavenumber is defined by 


k z = k c sin d = 



(VII. 7) 


where 0 is the angle of propagation relative to horizontal and k 0 is the reference 
wavenumber. 
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For the rest of this chapter, we assume that we have a point source, which is 
modeled in such way that [Ref. 35] 


= a(k 0 )S(f), 


(VII.8) 


where 


a(k 0 ) 


I iRp 

'I'Kko 


(YII.9) 


and S(f) is the amplitude spectrum of the pulse in the frequency domain. The MMPE 
uses a Hanning window in order to create the aforementioned spectrum, defined by 

c °s 2 ( mv ^ > \f~fc\< s f 


S(f ) = 


0 , 


BW 


1/- /cl > 2 


(VII.10) 


where f c is the center frequency of the transmitted band and BW is the transmission 
bandwidth. 

Now that there is motion of the source, Eq. (VII.6) should be stated in 
accordance with the transmission parameters (see [Ref. 35]) as 


i> (r = 0, fe, /) = h) - M, 


(VII.ll) 


where 


k z ,r = 


k 0 sin 6 

1 + ^cos‘ 


(VII. 12) 


The frequency in the enviromental frame is defined by the transmitted frequency, fc, 


with the assistance of Eq. (VII. 1), which determines the relationship between fc and 
fs{Q). Finally the source spectrum S(f) will take the form [Ref. 35] 


S(fc) = S 


f 


1 + ^ cos (6 - fc) 


(VII. 13) 


Finally the factor a has to be replaced by the corresponding transmitted value ax- 
One of the differences between the two cases, with and without Doppler, is the 
different size of the observed frequency bandwidth. In order to be able to compare 
results, the computational bandwidth of the non-Doppler case has to be set to the 
same size as the Doppler case. This is achieved by applying a weighting of zero to 
those frequency bins that are outside of the actual non-Doppler bandwidth. 
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B. MATLAB IMPLEMENTATION AND SIMULATION 
RESULTS 

For the tracking problem, the MATLAB implementation is almost the same 
as the implementation for the localization problem. However, there are two main 
differences. The first is the different manipulation of the MMPE data, since the code 
has now been changed in order to give us the effect of the Doppler. The second is an 
addition to the previous code, which has as inputs the coordinates of the buoys, the 
estimated coordinates of the source, the received signals correspon ding to each buoy, 
and the evaluated group speed (see Chapter 6, section 2). The computed outputs will 
then be the course and the speed of the target. Regarding the new mani p ula tion 
of the MMPE data, most of the procedure was shown in the previous section. So 
after re-mapping the original “MMPE” signal with the desired frequency resolution 
to account for the shift of the center frequency, we take the product of the real 
part of the re-mapped signal with the phasor exp( _i27r - fcTime ), where f c is the center 
frequency of the received signal, and Time is the re-mapped time vector with the 
desired resolution. Finally there is an interpolation of the pressure vector from its 
original time sp annin g vector Time to a fixed time spanning vector T im p fa which 
is the same for all received signals. 

For the computation of the course and the speed of the target, the procedure 
that we implemented was the following. First, we find the spectrum of each received 
signal and from the spectrum we extract the center frequency. Afterwards, with the 
coordinates of the buoy locations and the estimated target location, we extract the 
angles between the relative position vectors for each pair of the buoys. Then we find 
the relative Doppler frequency that corresponds to each pair of center frequencies of 
the received signals. Finally we estimate the course of the target for each buoy by the 
following equations. Note that for the simplicity of the problem we used only three 
receivers (N=3), otherwise the system of equations becomes very complex. 

Course( 1) = acos(A(l)) + B(l) (VII. 14) 
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Course{ 2 ) = acos(A( 2)) + B( 2) (VII.15) 

Course{ 3) = acos(A(3)) + .6(3) (VII.16) 


where 





A(l) 

(k%+ki) - jq 

2K,K 3 

(VII. 17) 


A( 2) 

(■Kf + *?) - Ki 

2K,K\ 

(VII. 18) 


A(3) 

(Ki + Ki) - KI 

2K X K 2 

(VII. 19) 

and 





B( 1) 

= 

(VII.20) 


B( 2) 

= 2 (^( 2 )+^( 3 )) 

(VII.21) 


B( 3) 

= i(*(3) + «l)) 

(VII.22) 

and 





Kr = 

A/ 12 ci 

sin (|(0(1) - 0(2))) 

(VII.23) 


k 2 = 

A/ 23 c| 

sin (|(0(2) - 0(3))) 

(VII.24) 


K z = 

A/ 31 C 5 

sin (^(0(3) - 0(1))) 

(VII.25) 


In Eqs. (VII.20) - (VII.22) and (VII.23) - (VII.25), the quantities 0(1), 0(2), 0(3) 
represent the angles in the horizontal plane of the position vectors of the buoys and in 
Eqs. (VII.23) - (VII.25) the variables A/i 2 , A/ 23 , A/ 3 i stand for the relative Doppler 
frequencies between the center frequencies of the received signals. In the end, the 
course of the target will be the average of the outcomes of Eqs. (VII. 14) - (VII. 16). 

Regarding the speed, it will be the average of the following calculations 


Speed( 1 ) 


_ A/12C _ 

2 / c sin ( ^(V-'K 2 ) ) sin (Course(l) — 


(VII.26) 
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Speed(2 ) 
Speed( 3) 


_ A/ 23 c _ 

2/c sin ( M-^ 3 ) ) sin (Course(2) - ) 

_ A/31C _ 

2/ c sin(- ^ 3 )-^ 1 ) ) sin (Course(Z) — < ^( 3 )+ < ^( 1 ) j 


(VII.27) 

(VTI.28) 


where / c is the theoretical center frequency of the transmitted signal. Since this 
quantity is unknown in our case, we substitute it with the average of the center 
frequencies of the received signals. This substitution will not significantly distort 
the result, since the shift of the center frequency due to Doppler is not large enough 
(low speeds) and the difference between the f c (transmitted) and the average of f c 
(received) is small. Figure 46 gives a geometrical overview of the problem. 



Figure 46. Top view of the Tracking problem. 

Some simulation results will follow, one from each enviro nm ent. For each 
case there will be figures indicating the geometric scenario and the solution that each 
method gives for the localization problem and also tables with Target’s Position Error 
in (m) and the Target’s motion Data, Course( 0 - 359 °) and Speed{m/s). 

1. Flat Botto m 

The characteristics of this range-independent environment remain the same as 
were presented in Chapters 4, 5 and 6 . In this case, the source is at depth 100 m and 
the 3 buoys are located at ranges 2, 3 and 5 km and at depth 30m as shown in Fig. 
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47. The source course and speed are 40° and 15m/s, respectively. Table XVI presents 
the results for the localization and tracking accuracy. 


Method 

Target’s Position Error in (m) 

MUSIC 

130.8 

Modified-MUSIC 

120.2 | 

MIN-NORM 

132.0 | 

PCLP 

131.1 

ESPRIT 

137.4 

Root-MUSIC 

137.1 

SCOT 

~ 1010.0 

PHAT 

1010.0 i 

Roth 

1293.0 i 

X-Corr 

53.6 

Target’s Motion 

Estimated 

Actual 

Course{ 0 — 359°) 

37.73 

40.00 

Speed (m/s) 

13.95 

15.00 


Table XVI. Target’s Motion & position Error for Flat bottom isospeed case (Number 
of receivers = 3). 
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2. Sound Channel 

As before, this range-independent environment retains the same properties as 
shown in Chapters 4, 5 and 6. For this case, the source is at depth 150m and the 3 
buoys are located at ranges 2 and 4km and at depth 40m as shown in Fig. 48. The 
source course and speed are 75° and lOm/s, respectively. Table XVII presents the 
results for the localization and tracking accuracy. 


Method 

Target’s Position Error in (m) 

MUSIC 

126.0 

Modified-MUSIC 

131.7 

MIN-NORM 

127.6 

PCLP 

126.3 

ESPRIT 

130.6 

Root-MUSIC 

129.4 

SCOT 

2234.5 1 

PHAT 

2234.5 

Roth 

4413.0 

X-Corr 

229.3 

Target’s Motion 

Estimated 

Actual 

Course (0 - 359°) 

75.17 

75.00 

Speed (m/s) 

11.41 

10.00 


Table XVII. Target’s Motion & position Error for Sound Channel case (Number of 
receivers = 3). 
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3. Shelf Break 

The properties of the current range-dependent environment are the same as 
described in previous Chapters. In this case, the source is at depth 250m and the 3 
buoys are located at ranges 3 and 4km and at depth 100m as shown in Fig. 49. The 
source course and speed are 100° and 8m/s, respectively. Table XVIII presents the 
results for the localization and tracking accuracy. 


Method 

Target’s Position Error in (m) 

MUSIC 

259.2 

Modified-MUSIC 

289.3 

MIN-NORM 

282.8 

PCLP 

275.6 

ESPRIT 

335.7 

Root-MUSIC 

303.5 

SCOT 

671.4 ~~ 1 

PHAT 

671.4 

Roth 

1947.6 

X-Corr 

185.2 

Target’s Motion 

Estimated 

Actual 

Course(0 - 359°) 

105.64 

100.00 

Speed (m/s) 

8.43 

8.00 


Table XVIII. Target’s Motion & position Error for Shelf Break case (Number of 
receivers = 3). 
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Figure 49. Tracking problem for Shelf Break case with N=3 sonobuoys. (a) Subspace 


methods, (b) Classical methods. 
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4. Internal Waves 

Finally this range-dependent environment retains the properties that we have 
already seen in our previous discussion. During the present simulation, the source is 
at depth 150m and the positions for the 3 buoys are at range 2km and at depth 50m 
as shown in Fig. 50. The source course and speed are 230° and 5m/s, respectively. 
Table XIX presents the results for the localization and tracking accuracy. 


Method 

Target’s Position Error in (m) 

MUSIC 

176.5 

Modifhed-MUSIC 

186.3 

MIN-NORM 

176.7 

PCLP 

176.2 

ESPRIT 

181.4 

Root-MUSIC 

181.3 

SCOT 

938.5 

PHAT 

938.5 

Roth 

1504.8 l 

X-Corr 

103.2 

Target’s Motion 

Estimated 

Actual 

Course (0 — 359°) 

228.77 

230.00 

Speed (m/s) 

4.21 

5.00 


Table XIX. Target’s Motion & position Error for Internal Waves case (Number of 
receivers = 3). 
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Figure 50. Tracking problem for Internal Waves case with N=3 sonobuoys. (a) 
Subspace methods, (b) Classical methods. 











Prom the presented results, we observe that the accuracy of the localization of 
the target has been decreased, compared with the cases where there was no movement 
of the source (zero Doppler). Especially for the classical methods, the error in posi¬ 
tion’s target was large enough to say that they failed to estimate the location of the 
source. One possible reason for the poor performance of the methods on the target 
localization is the approximation that we made during the manipulation of the data 
from MMPE when we were trying to set all signals in a fixed time spanning vector. 
On the other hand, the estimation of the target course and speed was quite accurate, 
or at least inside the acceptable limi ts of error, ±5° for the course and ±2m/s for the 
speed. 
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VIII. 


CONCLUSIONS AND FUTURE WORK 


A. CONCLUSIONS 

In this thesis, new techniques were developed to estimate the Time Difference 
of Arrival (TDOA) between two received signals at two seperate locations which 
contain the same transient but different noise. These techniques are based on subspace 
methods, a class of techniques based on the concept of signal and noise subspaces 
associated with the correlation matrix for a random process, and whose principal 
application is to locate narrowband hues in a spectrum or to estimate bearing of 
narrowband sources using array processing. As far as we know, the use of subspace 
methods to estimate TDOA for broadband transient sources is new. Our objectives 
were to implement the subspace methods in our problem, and to test them thoroughly 
comparing their performance with traditional methods of TDOA based on generalized 
cross-correlation (“classical methods”). The testing was conducted using synthetic 
data generated in MATLAB and data from an acoustic propagation model (MMPE). 
Tests were conducted under different scenarios using various transients, noise, or 
environmental conditions. 

The thesis went further to apply these methods in the localization of a target 
using the TDOA ranging algorithm and data from MMPE generated for various ocean 
environments. In this way, we tried to simulate realistic conditions in the ocean as 
much as possible. This second step was further expanded by the use of Doppler data 
from the MMPE algorithm to estimate not only the position but also course and 
speed; in other words, to track the source. 

From the various types of testing, we can conclude that the subspace methods 
provide good results in almost all cases. As a group, the subspace methods were 
more consistent in the accuracy of the results than the classical methods. This was 
observed not only in the TDOA testing but also in the localization and tracking 
experiments, where the subspace methods were more consistently accurate than the 
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classical methods. 

The plots produced by the subspace methods provide a clear and more easily 
identified correct peak than those corresponding to the classical methods and were 
not influenced as much by the noise. In addition some of these subspace methods 
can also be used to produce direct numerical estimates for the time delay without the 
need to search for a peak (ESPRIT - root MUSIC). 

Almost all methods (classical and subspace) performed satisfactorily in the 
localization problem. However accuracy tended to be more degraded for the range- 
dependent problems. This is not too surprising since the basic localization algorithm 
makes the assumption of a constant speed of sound. In the tracking experiments, 
the subspace methods had limited performance for estimation of target position but 
much better performance for estimation of course and speed. 

B. SUGGESTIONS FOR FUTURE DEVELOPMENT 

The results of this thesis show promising results for the subspace methods, 
and their application to target localization and tracking via TDOA estimation. The 
testing should be continued with more complex environments using acoustic models 
and also with real data from ocean experiments. A comparison of the subspace 
methods with other families of methods for TDOA (e.g., bicorrelation and bispectra 
[Ref. 6, 7, 8, 9]), or even a combination of these, would give a wider view of the 
possible solutions to the problem. Finally, implementation of a three-dimensional 
localization algorithm may give more realistic results and could assist in the tracking 
problem, provided that it can overcome obstacles such as multipath propagation. 



APPENDIX A. SYNTHETIC DATA 


This appendix presents sample figures and tables with summary results of the 
time delay estimation obtained from various forms of synthetic data. The transients 
that were used are the following: 

• Exponential 

• Sinusoidal 

• Modified-Exponential 

• Chirp 

• Modified-Chirp 

The definitions of the transients are given in Chapter 4 Eqs. (IV.2) - (IV.8). 
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Method 


MUSIC 


Modified MUSIC 


MIN-NORM 


ESPRIT 


Root-MUSIC 


SCOT 


PHAT 


Roth 


X-Corr 


Parameters 


Cov-Mat: 10 


Cov-Mat: 20 


Cov-Mat: 10 


Cov-Mat: 20 


Cov-Mat: 10 
Cov-Mat: 20 


Cov-Mat: 10 


Cov-Mat: 10 


Cov-Mat: 20 


Cov-Mat: 10 


Cov-Mat: 20 


#-Segs: 1 


#-Segs: 2 


#-Segs: 1 
#-Segs: 2 


#-Segs: 1 
#-Segs: 2 


#-Segs: 1 


#-Segs: 2 


SNR 15 dB 

SNR 10 dB 

SNR 5 dB 

101 

100 

106 j 

99 

99 

104 

101 

100 

107 

99 

99 

104 

100 

100 

110 

99 

98 

102 

100 

100 

109 

99 

98 

103 

100 

101 

107 

99 

98 

103 

101 

100 

106 

99 

99 

104 

102 

385 

5 

101 

102 

9 

102 

385 

5 

101 

202 

98 

-228 

315 

233 

101 

201 

9 

98 

107 

93 | 

101 

102 

109 


Table XX. Exponential transient: length L=20s; white noise (cr? = 2): actual 
TDOA=100s. 
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T(sec) T(sec) 

(a) 





T(sec) 


(b) 


Figure 51. Exponential transient: Actual Tdoa=100s; SNR=15dB; White-Noise (cr^). 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg¬ 
ments^. 
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(b) 

Figure 52. Exponential transient: Actual Tdoa=100s; SNR=15dB; White-Noise {a 2 0 ) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 
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Figure 53. Exponential transient: Actual Tdoa=100s; SNR=10dB; White-Noise (&%) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 

















































































































































































Figure 54. Exponential transient: Actual Tdoa=100s; SNR=10dB; White-Noise (a%) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 





























































































































Figure 55. Exponential transient: Actual Tdoa=100s; SNR=05dB; White-Noise (cr^) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 
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Figure 56. Exponential transient: Actual Tdoa=100s; SNR=05dB; White-Noise (o*) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, n umb er of seg 
ments=2. 































































































Method 

Parameters 

SNR 15 dB 

SNR 10 dB 

SNR 5 dB 

MUSIC 

Cov-Mat: 10 

-10 

-10 

-10 



-8 

-10 

Modified MUSIC 

Cov-Mat: 10 

-10 

-10 

-10 

Cov-Mat: 20 

-11 

-8 

-10 

MIN-NORM 

Cov-Mat: 10 

-10 

-10 

-9 

Cov-Mat: 20 

-10 

-9 

-10 

PCLP 

Cov-Mat: 10 

-10 

-10 

-9 

Cov-Mat: 20 

-10 

-9 

-10 

ESPRIT 

Cov-Mat: 10 

-10 

-10 

-9 

Cov-Mat: 20 

-10 

-9 

-14 

Root-MUSIC 

Cov-Mat: 10 

-10 

-10 

-9 

Cov-Mat: 20 

-10 

-8 

-10 

SCOT 

#-Segs: 1 

-10 

-10 

-76 

#-Segs: 2 

-10 

-11 


PHAT 


-10 

-10 

-76 

#-Segs: 2 

-10 

-15 

-91 

Roth 

#-Segs: 1 

-371 

-42 

-131 

#-Segs: 2 

-10 

-209 

169 

X-Corr 

#-Segs: 1 

-10 

-9 

-9 

#-Segs: 2 

-10 

-11 

-5 


Table XXI. Exponential transient: length L=50s; white noise = 2); actual 
TDOA=-10s. 
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T ( sec ) T(seo) 

(b) 

Figure 57. Exponential transient: Actual Tdoa=-10s; SNR=15dB; White-Noise ( c ? 0 ) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 

















Figure 58. Exponential transient: Actual Tdoa=-10s; SNR=15dB; White-Noise (a^) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 
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Figure 59. Exponential transient: Actual Tdoa=-10s; SNR=10dB; White-Noise (cr%) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, n umb er of seg 
ments=l. 

























































































































Figure 60. Exponential transient: Actual Tdoa=-10s; SNR=10dB; White-Noise (a^) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 
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Figure 61. Exponential transient: Actual Tdoa^ 
(a) Subspace methods, Covariance Size=10. (1 
ments=l. 


-10s; SNR=05dB; White-Noise (o*) 
) Classical methods, number of seg 






































































































































Figure 62. Exponential transient: Actual Tdoa=-10s; SNR=05dB; White-Noise (a%) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 
























































































































Method 

Parameters 



Damped 

Sin 

Chirp 

Damped 

Chirp 

MUSIC 

Cov-Mat: 10 

50 

49 

50 

50 

50 

Cov-Mat: 20 

50 

51 

50 

50 

50 

Modified-MUSIC 

Cov-Mat: 10 

50 

50 

50 

50 

50 

Cov-Mat: 20 

50 

50 

50 

50 

50 

MIN-NORM 

Cov-Mat: 10 

50 

49 

50 

50 

50 

Cov-Mat: 20 

50 

51 



50 

PCLP 

Cov-Mat: 10 

50 

49 

50 

50 

50 

Cov-Mat: 20 

50 

51 



50 

ESPRIT 

Cov-Mat: 10 

49 

49 

50 

50 

50 

Cov-Mat: 20 

50 




50 

Root-MUSIC 

Cov-Mat: 10 

50 

49 

50 

50 

50 

Cov-Mat: 20 

50 

51 

50 

50 

50 

SCOT 

#-Segs: 1 

50 

50 

50 

50 

50 

#-Segs: 2 

50 

-74 


50 

50 

PHAT 

#-Segs: 1 

50 

50 

50 

50 

50 

#-Segs: 2 

50 

-74 

50 

50 

50 

Roth 

#-Segs: 1 

-299 

-61 

50 

50 

50 

#-Segs: 2 

50 

— 

50 

50 

50 

X-Corr 

#-Segs: 1 

50 

50 

50 

50 

50 

#-Segs: 2 

50 

48 

50 

50 

50 


Table XXII. Transients: same length; white noise (o 2 = 2); SNR=15dB; actual 
TDOA=50s. 
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Figure 63. Sinusoidal transient: Actual Tdoa=50s; SNR=15dB; White-Noise (cr%) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 


























































































Figure 64. Sinusoidal transient: Actual Tdoa=50s; SNR=15dB; White-Noise (<%) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 






























































Figure 65. Sinusoidal transient: Actual Tdoa=50s; SNR=10dB; White-Noise (<r^) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 
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Figure 67. Sinusoidal transient: Actual Tdoa=50s; SNR=05dB; White-Noise (<t^) 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of seg 
ments=l. 




























































































































































Figure 68. Sinusoidal transient: Actual Tdoa=50s; SNR=05dB; White-Noise (a 2 0 ) 
(a) Subspace methods, Covariance Size=20. (b) Classical methods, number of seg 
ments=2. 













































































Method 

Parameters 

Exp 



Chirp 

Damped 

Chirp 

MUSIC 

Cov-Mat: 10 

51 

51 

50 

i 53 

49 

Cov-Mat: 20 

49 

51 

50 

50 

50 

Modified-MUSIC 

Cov-Mat: 10 

51 

52 

50 

53 

49 



51 

50 

50 

50 

MIN-NORM 



52 

50 

54 

49 

Cov-Mat: 20 

49 

51 

50 

50 

50 

PCLP 

itBgy/fliwni 


52 

50 

54 

49 



51 

50 

50 

50 

ESPRIT 

Cov-Mat: 10 

51 

52 

50 

54 

49 



51 

50 

50 

50 

Root-MUSIC 



52 

50 ! 

54 

49 

Cov-Mat: 20 

49 

51 

50 

50 

50 

SCOT 

#-Segs: 1 . 

50 

50 

50 

50 

50 

#-Segs: 2 

52 

-255 

50 

50 


PH AT 

^BDI 

50 

50 

50 

50 

50 

#-Segs: 2 

52 

48 

50 

50 

50 

Roth 

#-Segs: 1 

-502 

50 

108 

50 

50 | 


52 

-255 

50 

50 

50 

X-Corr 

#-Segs: 1 

50 

50 

50 

50 

50 

#-Segs: 2 

50 

48 

50 

50 

50 j 


Table XXIII. Transients: same length; white noise (<7„ = 2); SNR=10dB; actual 
TDOA=50s. 
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Figure 69. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=15dB; White- 
Noise (cr^). (a) Subspace methods, Covariance Size=10. (b) Classical methods, num¬ 
ber of segments=l. 














Figure 70. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=15dB; White- 
Noise (cr„). (a) Subspace methods, Covariance Size=20. (b) Classical methods, num¬ 
ber of segments=2. 



















Figure 71. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=10dB; White- 
Noise (o%). (a) Subspace methods, Covariance Size=10. (b) Classical methods, num¬ 
ber of segments=l. 

















































Figure 72. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=10dB; White- 
Noise (crl). (a) Subspace methods, Covariance Size=20. (b) Classical methods, num¬ 
ber of segments=2. 



























Figure 73. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=05dB; White- 
Noise (aj). (a) Subspace methods, Covariance Size=10. (b) Classical methods, num¬ 
ber of segments=l. 
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Figure 74. Damped Sinusoidal transient: Actual Tdoa=50s; SNR=05dB; White- 
Noise (al). (a) Subspace methods, Covariance Size=20. (b) Classical methods, num¬ 
ber of segments=2. 
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Method 

Parameters 

Exp 

Sin 

Damped 

Sin 

Chirp 


MUSIC 

Cov-Mat: 10 

51 

50 

49 

53 

51 


49 

52 

51 

48 

51 

Modified-MUSIC 

Cov-Mat: 10 

52 


49 

53 

51 

Cov-Mat: 20 

50 

51 

51 

48 

51 

MIN-NORM 

Cov-Mat: 10 

50 

50 

49 

57 

50 

Cov-Mat: 20 

49 

52 

51 

50 

50 

PCLP 

Cov-Mat: 10 

52 

50 

49 

56 

50 


49 

52 

51 

50 

50 

ESPRIT 

Cov-Mat: 10 

50 

50 

49 

53 

49 

Cov-Mat: 20 

50 

52 

51 

HI 

51 

Root-MUSIC 

Cov-Mat: 10 

49 

50 

49 

52 

50 

Cov-Mat: 20 

49 

52 

51 

49 

51 

SCOT 

#-Segs: 1 

92 

40 

52 

50 

50 

#-Segs: 2 

47 

mvm 

48 

50 

49 

PHAT 

#-Segs: 1 

92 

40 

52 

50 

50 

#-Segs: 2 

47 

46 

56 

50 

49 

Roth 

#-Segs: 1 

248 

478 

-204 

50 

50 

#-Segs: 2 



-216 

50 

49 | 

X-Corr 

#-Segs: 1 

53 

48 

46 

50 

50 

#-Segs: 2 

47 

46 

48 

50 

50_| 


Table XXIV. Transients: same length; white noise (o* = 2); SNR=05dB; actual 
TDOA=50s. 
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Figure 75. Chirp transient: Actual Tdoa—50s; SNR=15dB; White-Noise (cr^). (a) 
Subspace methods, Covariance Size=10. (b) Classical methods, number of seg¬ 
ments^. 
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Figure 76. Chirp transient: Actual Tdoa=50s; SNR=15dB; White-Noise (<%). (a) 
Subspace methods, Covariance Size=20. (b) Classical methods, number of seg- 
ments=2. 












Figure 77. Chirp transient: Actual Tdoa=50s; SNR=10dB; White-Noise (cr%). (a) 
Subspace methods, Covariance Size=10. (b) Classical methods, number of seg¬ 
ments^. 
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Figure 78. Chirp transient: Actual Tdoa=50s; SNR=10dB; White-Noise {a 2 0 ). (a) 
Subspace methods, Covariance Size=20. (b) Classical methods, number of seg- 
ments=2. 
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Figure 79. Chirp transient: Actual Tdoa=50s; SNR=05dB; White-Noise (<?%). (a) 
Subspace methods, Covariance Size=10. (b) Classical methods, number of seg¬ 
ments^. 
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Figure 80. Chirp transient: Actual Tdoa=50s; SNR=05dB; White-Noise (a%). (a) 
Subspace methods, Covariance Size=20. (b) Classical methods, number of seg- 
ments=2. 












Figure 81. Damped Chirp transient: Actual Tdoa=50s; SNR=15dB; White-Noise 
(<rj). (a) Subspace methods, Covariance Size=10. (b) Classical methods, number of 
segments=l. 
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Figure 82. Damped Chirp transient: Actual Tdoa=50s; SNR=15dB; White-Noise 
(oo)- ( a ) Subspace methods, Covariance Size=20. (b) Classical methods, number of 
segments=2. 
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Figure 83. Damped Chirp transient: Actual Tdoa=50s; SNR=10dB; White-Noise 
(a) Subspace methods, Covariance Size=10. (b) Classical methods, number of 
segments=l. 
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Figure 84. Damped Chirp transient: Actual Tdoa=50s; SNR=10dB; White-Noise 
(<%)• ( a ) Subspace methods, Covariance Size=20. (b) Classical methods, number of 
segments=2. 
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Figure 85. Damped Chirp transient: Actual Tdoa=50s; SNR=05dB; White-Noise 
(aj). (a) Subspace methods, Covariance Size=10. (b) Classical methods, number of 
segments=l. 
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Figure 86. Damped Chirp transient: Actual Tdoa=50s; SNR=05dB; White-Noise 
(O- (a) Subspace methods, Covariance Size=20. (b) Classical methods, number of 
segments=2. 
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